Open Access
Add to BibSonomy
Abstract
We propose a method to design double smooth freeform surfaces applied in beam shaping with a ray mapping method in the paper. We couple the calculation of ray mapping and the construction of freeform surfaces to approach the surface normal field integrability condition based on the symplectic flow mapping scheme. In this paper, the incident beam wavefront is not limited to be planar or spherical. Several challenging design examples are presented that include transforming a circular Gaussian beam to an unconventional beam with variously shaped contour, and transforming an elliptic beam to a convergent beam with complex irradiance distribution in non-paraxial regime. The results show the high efficiency and feasibility of the proposed method in designing freeform optics for beam shaping applications.
© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
Beam shaping is the process of redistributing the irradiance and phase of a beam of optical radiation, which is playing a critical role in many laser applications such as optical data storage, laser printing, lithography and materials processing [1,2]. The problem can be well addressed by using lenses or mirrors with two elaborately designed surfaces when geometrical optics is applicable [3]. The advantages of using refractive and reflective optics in beam shaping include high energy transfer efficiency, compact structure and easy implementation.
In earlier researches, lots of studies of refractive or reflective beam-shaping system design for controlling irradiance and wavefront is limited to rotationally symmetric optics [4–9]. For more general cases, designing optical elements without assumptions of any symmetry is quite more challenging. One way to achieve the freeform surface design is converting the problem to an elliptic nonlinear partial differential equation (PDE) [10,11] and then solving it numerically. However, obtaining a stable solution of this nonlinear equation for different design requirements is still a theoretically and computationally complex task. An alternative approach to calculate two freeform surfaces for irradiance and phase control is based on the mathematically rigorous geometric supporting quadric method (SQM) [12] which is very efficient and robust for the design of freeform illumination optics proposed by Oliker. Oliker also has done many other beneficial works of rigorous mathematical solutions for freeform surface design controlling irradiance and phase of a collimated beam [13–15].
A simple and flexible approach to address the problem is based on the ray mapping method. This method is generally implemented in two steps by firstly establishing a mapping between the incident and output beams and then constructing the freeform surfaces according to the calculated ray mapping. In this method, the energy transportation from incident and output beams is usually regarded as optimal mass transportation (OMT) problem with a quadric cost function [16–18]. Nevertheless, the ray mapping method based on OMT with quadric cost function can satisfy the surface normal integrability condition only in paraxial approximation [18–20].
In [21] and [22], the authors proposed a method based on OMT with non-quadric cost function and converted the calculation of ray mapping to linear assignment problems. This approach calculates the ray mapping for collimated beams in non-paraxial cases and then reconstructs the freeform surfaces with B-spline surface and SQM method. The method is very efficient for achieving a piecewise smooth freeform surface for domains with non-smooth boundaries and disconnected regions.
Current development of designing freeform surfaces for irradiance and phase control particularly focuses on collimated incident beams. For beams with arbitrary wavefront, designing double freeform surfaces for beam shaping becomes quite a trouble for most widely published approaches. Another issue arises in the ray mapping method. The traditional ray mapping method generally implements the calculation of ray mapping and construction of freeform surface independently. This feature can bring difficulties to achieve the ray mapping that yields an integrable normal field because the normal field can’t be evaluated accurately without surface information. Our paper is based on the symplectic flow ray mapping scheme proposed by Karel Desnijder et al. [23]. This method calculates the ray mapping by directly eliminating the curl of surface normal for arbitrary wavefront input beams which may have the potential to provide a general methodology for designing freeform beam shaping systems. It is worth to be noted that the method described in [23] is applied to far-field approximation, for which the surface normal field can be directly calculated from the ray mapping. For beam shaping and the cases that the size of freeform lens cannot be neglected, it has to start the algorithm using far-field approximation and repeat the algorithm many times, which is cumbersome and increases computing time. In this paper, we couple the calculation of ray mapping and construction of freeform surface in the same iteration step, which can tackle the calculation of integrable ray mapping of beam shaping effectively.
There are several researches that provide similar objectives or thoughts of this paper to achieve freeform illumination system design. In [24], the authors provide a ray mapping method for designing double freeform surfaces that can generate prescribed irradiance at two target planes from a divergent beam. In [25], the authors combine the construction of freeform surface with optimal transport ray mapping to design freeform reflector for off-axis illumination. Feng recently proposed an iterative wavefront tailoring (IWT) method [26] which can solve complex problems by solving a simplified Monge-Ampére equation and operating iteration steps for a point-like source. The method [26] couples the surface construction in every iteration step to revise the output wavefront and approach integrable ray mapping, which also has the potential to provide a universal methodology to tackle beam shaping problem.
There are two main purposes of our work: (1) this paper aims to provide a method to design double smooth freeform surfaces for controlling the irradiance and phase of beams with no limitation of planar or spherical input wavefront; (2) this method is proposed to efficiently achieve the ray mapping that yields integrable surface normal field for the cases that the size of freeform lens cannot be neglected. The paper is organized as follows. In Section 2, we give a formulation of the freeform beam shaping design problem. In Section 3, we demonstrate the method detailedly in three subsections including the initial mapping calculation, surface construction and iteration process. In Section 4, several challenging design examples of transforming incident Gaussian beams to unconventional beams are presented to prove the feasibility of the proposed method for freeform optics design in beam shaping applications. In Section 5, a brief conclusion is presented.
2. Formulation of the problem
We focus on the freeform refractive system for beam shaping in this paper. For reflective system design, the method is exactly the same as refractive system design approach demonstrated in this paper. The geometrical structure of the freeform surface design for beam shaping is presented schematically in Fig. 1. The source and target plane are located in a medium of refractive index $n_o$, and the refractive index of the medium between the double freeform surfaces is $n_i$. If $n_o = 1$ and $n_i\,>\,1$, Fig. 1 represents a single lens with two freeform surfaces, and for the case of $n_i = 1$ and $n_o\,>\,1$ the layout represents a two freeform lens system. Let $(x_1,x_2,z) \in \boldsymbol {R}^3$ denote the Cartesian coordinates with $z$ the horizontal coordinate and $\boldsymbol {x} = (x_1,x_2) \in \boldsymbol {R}^2$ represent the coordinates of source plane located at $z$ = 0. The coordinates of target plane situated at $z = l$ are defined as $\boldsymbol {y} = (y_1,y_2) \in \boldsymbol {R}^2$. $\mathcal {X}$ and $\mathcal {Y}$ denote the domains of source and target irradiance.
The purpose is to achieve the prescribed irradiance $g(\boldsymbol {y})$ at the target plane and desired output wavefront corresponding to the incident beam with irradiance $f(\boldsymbol {x})$ at source plane and input wavefront $\boldsymbol {W}$ by using two elaborately designed freeform surfaces. We assume that both freeform surfaces are ideal surfaces with no energy losses in the refraction.
3. Method
3.1 Calculation of initial mapping
The first step is to establish an initial mapping that satisfies the energy conservation between source and target plane. Using the method proposed in [27], we create a virtual rectangular starting grid and calculate the mapping from the starting grid to source and target separately. Assume that $(\xi ,\eta ) \in \boldsymbol {R}^2$ denote the coordinates of the virtual starting grid and $\mathcal {H}$ denotes the rectangular domain of the virtual uniform irradiance. The irradiance distribution of the starting grid is defined as $h(\xi ,\eta )$. The optimal mass transport mapping can be a favorable starting point because it is somewhat close to the final solution and also has been studied well in the widely published paper [16–18]. Let us consider the calculation of mapping between the virtual starting grid and the target. Let $\boldsymbol {m}$ represents the mapping from starting grid $(\xi ,\eta )$ to target $(y_1,y_2)$, which means $\boldsymbol {y}=\boldsymbol {m}(\xi ,\eta )$. The computation process on the source is the same as the mapping calculation for the target. The problem can be converted to the following equations:
where $D\boldsymbol {m}$ denotes the Jacobi matrix of mapping $\boldsymbol {m}$. Eq. (1) is the energy conversation equation and Eq. (2) represents that the boundary of the starting grid is mapped to the boundary of the target. For optimal transportation problem with a quadric cost function, the mapping is the unique gradient of a convex function, which means $\boldsymbol {m}=\nabla {u}$. Substituting $\boldsymbol {m}=\nabla {u}$ in Eq. (1), the equation is transformed into the standard Monge-Ampére equation. Several numerical algorithms can be applied to solve this standard Monge-Ampére equation such as Newton’s method [28–30] or converting the problem into an evolution problem relying on the gradient of time [16,31,32].In this section, we employ a least-squares method proposed by Prins et al. [33] which is very efficient and robust for the optimal transport problem. This method has also been successfully used in freeform illumination design with OMT ray mapping [34]. The method uses the fact that the mapping solving the optimal transport problem is equal to the gradient of the convex solution of Monge-Ampére equation with boundary condition Eq. (2). Therefore, the Jacobi matrix of $\boldsymbol {m}$ equals a real symmetric positive semidefinite matrix $\boldsymbol {P}(\xi ,\eta )$. Then, the method converts the problem to the minimization problems for the following three objective functions:
3.2 Double freeform surfaces construction through the point-by-point procedure
In this subsection, we give a brief overview of the point-by-point freeform surface construction method proposed by Feng et al. [16]. The method is chosen due to its speediness and efficiency in freeform surface construction during the iteration process. We will demonstrate the iteration method in the next subsection.
According to the calculated mapping and the prescribed wavefront, we can get two bundles of rays of incident and output beam. Let $I_{i,j}$ denote the incident ray sequence and $O_{i,j}$ denotes the output ray sequence, and $m,n$ denote the maximum subscript of the sequence. We need to define two fixed points on the first freeform surface and the second freeform surface through the first ray path of the ray bundles. Let $P_{1,1}$ and $Q_{1,1}$ represents the two fixed points on the two freeform surfaces respectively. The vector $\boldsymbol {P_{1,1}Q_{1,1}}$ denotes the first refracted ray between the double surfaces. We can calculate the optical path length (OPL) via the two fixed points. Assume that we have obtained two points $P_{i,1}$ and $Q_{i,1}$. We can firstly evaluate an approximate value of $P_{i+1,1}$ by intersecting the input ray $I_{i+1,1}$ with the tangent plane of $P_{i,1}$. Let $P'_{i+1,1}$ denotes the approximate point. Then, an approximate point on the second surface $Q'_{i+1,1}$ can be calculated by equaling the OPL. An approximate normal vector can be obtained according to $P'_{i+1,1}$ and $Q'_{i+1,1}$ which is denoted as $N'_{i+1,1}$. Next, we compute an intermediate point $P_{i+1/2,1}$ by intersecting $I_{i+1/2,1}$ with the tangent plane of $P_{i,1}$, where $I_{i+1/2,1}$ can be obtained as $(I_{i+1,1}+I_{i,1})/2$. After that, the point $P_{i,1}$ can be computed by intersecting $I_{i+1,1}$ with the plane through $P_{i+1/2,1}$ and perpendicular to $N'_{i+1,1}$. Then, we can get $Q_{i+1,1}$ based on the constant OPL condition. The iteration step sketch diagram is presented in Fig. 3(a).
Implementing the above procedure from $i=1$ to $m$, we can get the first curve denoted by $j=1$ on the two freeform surfaces. In the same way, if we have obtained the points $P_{i,j}$ and $Q_{i,j}$, we can calculate the point $P_{i,j+1}$ and $Q_{i,j+1}$ using the above method. Every point on curve $j=1$ is regarded as the starting point to perform the iteration step from $j=1$ to $n$. Then, we can get all points on the freeform surface via the iteration process. The conceptual diagram is shown in Fig. 3(b). The detailed method for the freeform construction can refer to [16].
3.3 Achieving integrable mapping via symplectic transformation
In this subsection, we couple the construction of freeform surface into the iteration process of ray mapping calculation via a symplectic flow. Compared to the traditional ray mapping method, whereby the establishment of ray mapping and the formulation of freeform surface are independent of each other, this approach calculates the ray mapping that satisfies the surface normal integrability condition with constructing freeform surface at the same time. This feature allows the method to tackle the freeform optics design for beam shaping efficiently.
From the process demonstrated in Section 3.1, we can obtain an initial curl-free mapping from a rectangular starting grid to a target. Assume that $\boldsymbol {m}$ denotes this mapping. According to the initial mapping from source to target, we can construct the two freeform surfaces and complete the design. That is the freeform optics design method known as ray mapping with optimal transport, which unfortunately can lead to a satisfactory result only in paraxial approximation. The aim is to calculate the ray mapping that lead to the integrable normal field of the first freeform surface, whereby we can construct the second freeform surface using the principle of equal optical path length (OPL). The first surface normal field can be calculated using Snell’s law:
where R and I are defined in Fig. 1. According to [23], to satisfy the normal integrability condition $\bf {N}\cdot (\nabla \times \bf {N})=0$, we only need to satisfy: where $\bf {N'}$ is defined as $\bf {N}/(\bf {N}\cdot \bf {I})$. Substituting Eq. (6) in Eq. (7), the left branch of Eq. (7) can be converted to:In order to find an integrable mapping, we need to alter the initial ray mapping in a way that Eq. (9) is satisfied. As demonstrated in [23], this can be done via a symplectic transformation of the coordinates $(\xi ,\eta )$. However, the refracted rays $\textbf {R}$ cannot be evaluated directly without surface information. In this paper, we couple the construction of freeform surface into the iterative process to find an integrable ray mapping for beam shaping. The point-by-point surface construction approach introduced in Section 3.2 is employed in the iterative process.
We defined the left branch of Eq. (9) as the cost function $C(\xi ,\eta )=(\nabla \times \textbf {R}')\cdot \textbf {I}$. By applying the formal definition of the curl of vector field, we write the integral in terms of partial of $\textbf {R}'$ and W:
The iteration process is presented in Fig. 4(a). For every iteration step, we need to construct the double freeform surfaces using the updated target mapping $\textbf {m}(\xi ,\eta )$ with Feng’s method. And we will get two sequences of points $P_{i,j}$ and $Q_{i,j}$ ($i=1,2,\ldots ,m;j=1,2,\ldots ,n$) of the two freeform surfaces as shown in Fig. 4(b). The vector $\boldsymbol {P_{i,j}Q_{i,j}}$ is the desired refracted ray vector. If the ray mapping can not yield an integrable normal vector field of the first freeform surface, the Eq. (9) is not satisfied and the value of cost function Eq. (10) is not zero. Therefore, the desired refracted rays will deviate from the actual refracted rays. The cost function and the Hamiltonian are calculated according to the desired refracted rays normalized vector field $\textbf {R}$. Then we can update the phase space coordinates $(\xi ,\eta )$ via Eq. (14). As the cost function decreases, the desired refracted rays and the surface normal get close to the actual refracted rays and surface normal. The process stops if the value of the cost function is close to zero. Note that the energy conservation holds at every point of phase space in coordinates $(\xi ,\eta )$ during the iteration step and that is a consequence of the famous Liouville’s theorem: phase-space volume is preserved. Notice that the Hamiltonian should be reframed in every iteration step according to the updated target mapping. When the calculation of the integrable mapping is completed, the two freeform surfaces are constructed at the same time. The convergence of the algorithm and the running time will be demonstrated in the next section with several challenging design examples.
4. Design examples
In this section, several design examples are presented to show the feasibility of the proposed method. We mainly focus on the two freeform lenses design problem for the case of $n_i = 1$ and $n_o\,>\,1$ labeled in Fig. 1. The refractive index of element material is 1.4856, which is the refractive index of fused silica at the wavelength of 308 nm.
4.1 Beam shaping with planar incident and output wavefront
Let us first consider the design of two plano-freeform lenses for collimated beam shaping as shown in Fig. 5, where $D$ represents the diameter of incident beams and $l$ denotes the distance between the plane surfaces of the first and second freeform lens. The incident beam with circular cross-section and Gaussian irradiance distribution $E(\boldsymbol {x})=e^{-2(x_1^2+x_2^2)/r_0^2}$,$\boldsymbol {x}\in \mathcal {X}=\{(x_1,x_2)|x_1^2+x_2^2 \leq (D/2)^2\}$ is considered. The parameters of the incident beam are determined as $r_0$=1 mm,$D$=2 mm.
The first example is considered as transforming the incident collimated Gaussian beam to a beam with rectangular cross-section of 4 mm$\times$3 mm and uniform irradiance. $l$ is determined as 7 mm. The $z$ coordinates of fixed point $P_{1,1}$ and $Q_{1,1}$ in Fig. 3 are set to be 1 mm and 6 mm. The method is implemented by setting 121$\times$121 virtual starting grid in Fig. 2. Applying the method proposed in Section 3, we obtain the final system presented in Fig. 6(a). Figure 6(a) shows that the irradiance distribution remains almost the same during the propagation of the beam. The comparison of the design from the OMT mapping and integrable mapping is specifically presented in Fig. 6(b) and 6(c). As shown in Fig. 6(c), the performance of irradiance distribution generated by lenses designed from integrable mapping is much better than that from OMT mapping in the non-paraxial regime.
Fig. 6. (a) The designed beam shaping system using integrable mapping; (b,c) the comparison between OMT mapping and integrable mapping.
Figure 7 presents the irradiance distributions in the planes $z$=10 mm and $z$=50 mm. The comparisons of the simulated results and prescribed irradiance distributions at different distances are also demonstrated in Fig. 7. The results show that the irradiance distribution hold up during the propagation of the output beam to a long distance. The iteration step is implemented 30000 times. We evaluate the value of the function $H$ in Eq. (12). The evaluation of merit function is determined as
where $m,n$ denote the number of starting virtual grids which are 121$\times$121 in this example. The convergence of the merit function is shown in Fig. 8. The merit function value is close to zero after 10000 times of iteration steps. The running time is also presented in Fig. 8. In this design example, we perform the iteration step Eq. (14) for 30000 times. The running time in Matlab is no more than 150 seconds on an Intel Core i7 8700k with 32GB RAM. The running time of surface construction accounts for 34$\%$ of one total iteration step.Fig. 7. (a) The irradiance distribution and row irradiance in the plane $z$=10 mm; (b) The irradiance distribution and row irradiance in the plane $z$=50 mm.
To demonstrate the efficiency of the proposed method, we have designed several challenging systems to transform the incident beam to unconventional beams with arbitrary contour and complex irradiance. Figure 9 presents the irradiance distributions generated by the designed systems for uniform beams with hexagonal and fan-shaped contour and beam with prescribed irradiance distribution ”OEI”. In these examples, $l$ is determined as 7 mm. The $z$ coordinates of fixed points $P_{1,1}$ and $Q_{1,1}$ in Fig. 3 are set to be 1 mm and 9 mm. All the configurations of the design examples are non-paraxial systems. As presented in Fig. 9, the irradiance distributions holds up during the propagation of the output beams. The results show that the designed systems achieve good irradiance and wavefront control for challenging design objectives. The time complexity of the iteration algorithm is $O(n^2)$ for an $n\times n$ grid. We implement the ”OEI” example on a finer grid 241$\times$241 and the running time for 20000 iteration steps is about 726s.
Fig. 9. Irradiance distributions of uniform hexagonal beam, fan-shaped beam and rectangular beam with prescribed irradiance ”OEI” generated by lenses in the plane $z$=12 mm and $z$=30 mm.
4.2 Beam shaping with non-planar incident and output wavefront
The next design example is redistributing a beam with a non-planar wavefront to a convergent beam with complex irradiance. The design sketch diagram is presented in Fig. 10(a). We predefine the entrance surface as an elliptic paraboloid $z=z_0(1-x^2/x_0^2-y^2/y_0^2)$ and let a collimated circular Gaussian beam incident on the entrance surface. The collimated Gaussian beam is the same as the first design example. The parameters of the entrance surface are determined as: $x_0$=2mm, $y_0$=1.5mm and $z_0$=2mm. Refracted by the entrance surface, the incident beam irradiance and wavefront are presented in Fig. 10(b) and 10(c).
Fig. 10. (a) The design sketch diagram of non-planar input wavefront beam; (b) the incident beam irradiance; (c) the incident wavefront.
The output beam is determined as a convergent beam whose center is located at $z$=30 mm. The target plane is at $z$=23 mm and the prescribed irradiance distribution is defined as a letter ”B” with a uniform rectangular background. The irradiance ratio is 1 (background) to 2.4 (letters). To determine the system parameters, we trace the first ray of the ray mapping denoted by $\textbf {r}_{1,1}$ traveling from the entrance surface to the last surface. The path is set to be 2 mm between the entrance surface to the first freeform surface along the first ray. The path of $\textbf {r}_{1,1}$ between the two freeform surfaces is 14 mm. This path also represents that the distance of fixed points $P_{1,1}$ and $Q_{1,1}$ in Fig. 3 is 14 mm. Then, we implement the method introduced in Section 3 to accomplish the design. As presented in Fig. 11(a), the ray mapping yields an integrable surface normal vector field between the source and target. The simulated irradiance distribution is shown in Fig. 11(b). The column and row irradiance as well as the comparisons between the simulated results and desired results are presented in Fig. 11(c).
Fig. 11. (a) The integrable mapping of source and target; (b) the simulated irradiance distribution at the target plane of the designed system; (c) the row and column irradiance distributions as well as the comparisons of the simulated results and the desired results.
The iteration step is implemented for 20000 times with 121$\times$121 starting grids. The value of the merit function is defined as Eq. (15). The convergence of the merit function is shown in Fig. 12. The merit function value is close to zero after 5000 times of iteration steps. The running time is also presented in Fig. 12. In this design example, we perform the iteration step Eq. (14) for 20000 times. The running time in Matlab is about 182 seconds on an Intel Core i7 8700k with 32GB RAM. The result of integrable ray mapping has been shown in Fig. 11(a).
Because of the convergent output wavefront, the relative irradiance distribution will remain the same during the propagation of the output beam. We simulated the irradiance distributions at different distances in the direction of the output beam propagation. The designed system and relative irradiance distributions are shown in Fig. 13. The results show that the designed systems achieve good irradiance and wavefront control in this challenging design problem.
In this method, the problem is reduced to find a mapping that minimizing the curl of normal field $\bf {N'}$ ($\bf {R'}$ in this paper) everywhere. And we can reduce the curl by twisting the vector field in the right direction. Since we start from a regular grid and implement the symplectic flow on coordinates $(\xi ,\eta )$ to ensure that the $(\xi ,\eta )$ grid remains uniform in the motion. Note that the Hamiltonian is recalculated in every iteration according to the updated target mapping. This process does not ensure that we can always approach the integrability condition, which means the cost function may not converge to zero all the time. However, we can always get much better results than the initial ray mapping and obtain a preferable freeform system without solving complex nonlinear partial differential equations.
5. Conclusion
In this paper, we have proposed a general method for designing double smooth freeform surfaces to control the irradiance and wavefront of incident beams with no limitation of input plane wavefront. We approach the surface normal vector integrability condition by directly controlling the refracted rays’ vector field using Eq. (9). Compared to the traditional ray mapping approach, the presented method couples the calculation of ray mapping and construction of freeform surface in the same process, which solves the integrable ray mapping for beam shaping efficiently. We have presented several challenging design examples of redistributing the irradiance and wavefront of incident Gaussian beams to complex irradiance distributions. The performance of the designed systems has demonstrated the efficiency and feasibility of the proposed method.
Funding
National Natural Science Foundation of China (61805088); Fundamental Research Funds for the Central Universities (2019kfyRCPY083, 2019kfyXKJC040).
Acknowledgment
The authors thank Karel Desnijder for valuable discussions of the numerical algorithm to calculate the ray mapping with surface normal field integrability condition satisfied.
References
1. F. M. Dickey and S. C. Holswade, Laser Beam Shaping: Theory and Techniques (Marcel Dekker, 2000).
2. F. M. Dickey, S. C. Holswade, T. E. Lizotte, and D. L. Shealy, Laser Beam Shaping Applications (Taylor Francis Group, 2006).
3. V. Galindo, “Design of dual-reflector antennas with arbitrary phase and amplitude distributions,” IEEE Trans. Antennas Propag. 12(4), 403–408 (1964). [CrossRef]
4. B. R. Frieden, “Lossless conversion of a plane laser wave to a plane wave of uniform irradiance,” Appl. Opt. 4(11), 1400–1403 (1965). [CrossRef]
5. J. L. Kreuzer, “Coherent light optical system yielding an output beam of desired intensity distribution at a desired equiphase surface,” (1969). US Patent 3,476,463.
6. P. W. Rhodes and D. L. Shealy, “Refractive optical systems for irradiance redistribution of collimated radiation: their design and analysis,” Appl. Opt. 19(20), 3545–3553 (1980). [CrossRef]
7. J. A. Hoffnagle and C. M. Jefferson, “Design and performance of a refractive optical system that converts a gaussian to a flattop beam,” Appl. Opt. 39(30), 5488–5499 (2000). [CrossRef]
8. H. Ma, Z. Liu, P. Jiang, X. Xu, and S. Du, “Improvement of galilean refractive beam shaping system for accurately generating near-diffraction-limited flattop beam with arbitrary beam size,” Opt. Express 19(14), 13105–13117 (2011). [CrossRef]
9. D. L. Shealy and S.-H. Chao, “Geometric optics-based design of laser beam shapers,” Opt. Eng. 42(11), 3123–3139 (2003). [CrossRef]
10. Y. Zhang, R. Wu, P. Liu, Z. Zheng, H. Li, and X. Liu, “Double freeform surfaces design for laser beam shaping with monge–ampère equation method,” Opt. Commun. 331, 297–305 (2014). [CrossRef]
11. 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-ampère type equation,” J. Opt. 18(12), 125602 (2016). [CrossRef]
12. V. Oliker, L. L. Doskolovich, and D. A. Bykov, “Beam shaping with a plano-freeform lens pair,” Opt. Express 26(15), 19406–19419 (2018). [CrossRef]
13. V. Oliker, “Optical design of freeform two-mirror beam-shaping systems,” J. Opt. Soc. Am. A 24(12), 3741–3752 (2007). [CrossRef]
14. V. Oliker, “On design of free-form refractive beam shapers, sensitivity to figure error, and convexity of lenses,” J. Opt. Soc. Am. A 25(12), 3067–3076 (2008). [CrossRef]
15. V. Oliker, “Designing freeform lenses for intensity and phase control of coherent light with help from geometry and mass transport,” Arch. Ration. Mech. Anal. 201(3), 1013–1045 (2011). [CrossRef]
16. 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]
17. Z. Feng, B. D. Froese, C.-Y. Huang, D. Ma, and R. Liang, “Creating unconventional geometric beams with large depth of field using double freeform-surface optics,” Appl. Opt. 54(20), 6277–6281 (2015). [CrossRef]
18. C. Bösel and H. Gross, “Ray mapping approach for the efficient design of continuous freeform surfaces,” Opt. Express 24(13), 14271–14282 (2016). [CrossRef]
19. C. Bösel, N. G. Worku, and H. Gross, “Ray-mapping approach in double freeform surface design for collimated beam shaping beyond the paraxial approximation,” Appl. Opt. 56(13), 3679–3688 (2017). [CrossRef]
20. L. L. Doskolovich, A. A. Mingazov, D. A. Bykov, E. S. Andreev, and E. A. Bezus, “Variational approach to calculation of light field eikonal function for illuminating a prescribed region,” Opt. Express 25(22), 26378–26392 (2017). [CrossRef]
21. L. L. Doskolovich, D. A. Bykov, E. S. Andreev, E. A. Bezus, and V. Oliker, “Designing double freeform surfaces for collimated beam shaping with optimal mass transportation and linear assignment problems,” Opt. Express 26(19), 24602–24613 (2018). [CrossRef]
22. D. A. Bykov, L. L. Doskolovich, A. A. Mingazov, E. A. Bezus, and N. L. Kazanskiy, “Linear assignment problem in the design of freeform refractive optical elements generating prescribed irradiance distributions,” Opt. Express 26(21), 27812–27825 (2018). [CrossRef]
23. K. Desnijder, P. Hanselaer, and Y. Meuret, “Ray mapping method for off-axis and non-paraxial freeform illumination lens design,” Opt. Lett. 44(4), 771–774 (2019). [CrossRef]
24. Z. Feng, B. D. Froese, R. Liang, D. Cheng, and Y. Wang, “Simplified freeform optics design for complicated laser beam shaping,” Appl. Opt. 56(33), 9308–9314 (2017). [CrossRef]
25. C. Gannon and R. Liang, “Ray mapping with surface information for freeform illumination design,” Opt. Express 25(8), 9426–9434 (2017). [CrossRef]
26. Z. Feng, D. Cheng, and Y. Wang, “Iterative wavefront tailoring to simplify freeform optical design for prescribed irradiance,” Opt. Lett. 44(9), 2274–2277 (2019). [CrossRef]
27. 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]
28. R. Wu, Y. Zhang, M. M. Sulman, Z. Zheng, P. Benítez, and J. C. Miñano, “Initial design with l 2 monge-kantorovich theory for the monge–ampère equation method in freeform surface illumination design,” Opt. Express 22(13), 16161–16177 (2014). [CrossRef]
29. B. D. Froese and A. M. Oberman, “Convergent filtered schemes for the monge–ampére partial differential equation,” SIAM J. on Numer. Analysis 51(1), 423–444 (2013). [CrossRef]
30. 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]
31. S. Haker, L. Zhu, A. Tannenbaum, and S. Angenent, “Optimal mass transport for registration and warping,” Int. J. computer vision 60(3), 225–240 (2004).
32. M. M. Sulman, J. Williams, and R. D. Russell, “An efficient approach for the numerical solution of the monge–ampère equation,” Appl. Numer. Math. 61(3), 298–307 (2011). [CrossRef]
33. C. Prins, R. Beltman, J. ten Thije Boonkkamp, W. IJzerman, and T. W. Tukker, “A least-squares method for optimal transport using the monge-ampere equation,” SIAM J. on Sci. Comput. 37(6), B937–B961 (2015). [CrossRef]
34. S. Wei, Z. Zhu, Z. Fan, Y. Yan, and D. Ma, “Multi-surface catadioptric freeform lens design for ultra-efficient off-axis road illumination,” Opt. Express 27(12), A779–A789 (2019). [CrossRef]
35. L. Casetti, “Efficient symplectic algorithms for numerical simulations of hamiltonian flows,” Phys. Scr. 51(1), 29–34 (1995). [CrossRef]
