## Abstract

We propose an improved double freeform-optical-surface design method for shaping a prescribed irradiance distribution whilst forming a desired wavefront from a given incident beam. This method generalizes our previous work [Opt. Exp. 21, 14728-14735 (2013)] to tackle non-separable beam irradiances. We firstly compute a proper ray mapping using an adaptive mesh method in the framework of the L^{2} Monge-Kantorovich mass transfer problem. Then, we construct the two freeform optical surfaces according to this mapping using a modified simultaneous point-by-point procedure which is aimed to minimize the surface errors. For the first surface, the modified procedure works by firstly approximating a value to the next point by only using the slope of the current point and then improving it by utilizing both slopes of the two points based on Snell’s law. Its corresponding point on the second surface can be computed using the constant optical path length condition. A design example of producing a challenging irradiance distribution and a non-ideal wavefront demonstrates the effectiveness of the method.

© 2013 Optical Society of America

## 1. Introduction

The problem we consider here is how to shape a prescribed irradiance distribution whilst forming a desired wavefront (or phase) from a given incident beam in the geometrical optics approximation. Generally, this involves designing at least two reflective or refractive freeform optical surfaces for increasing need of non-rotational symmetry transformations. The designs can play great roles in laser beam shaping applications and many other fields including illumination, solar energy concentration and reflector antennas.

Even the problem of only controlling the irradiance using one or two freeform optical surfaces is quite difficult. A variety of methods has been proposed to tackle this problem theoretically and numerically (see e.g. [1–13].). A second requirement for controlling the wavefront using additional freeform optical surfaces makes the design more demanding. Here, we only consider designing two freeform optical surfaces for the two requirements.

A direct determination of both freeform optical surfaces can lead to a very hard mathematical problem, since it needs to solve complicated simultaneous equations governing both the irradiance and wavefront [14]. Glimm and Oliker illustrated that a two-reflector system for reshaping an arbitrary irradiance whilst keeping the beam planar can also be solved as a variational problem within the framework of the L^{2} (quadratic cost) Monge-Kantorovich mass transfer problem [15]. This approach can also be extended to the design of refractive systems [16]. Rubinstein and Wolansky showed that a single lens with two freeform surfaces for shaping arbitrary collimated beams can also be found by solving a variational problem related to the weighted least action [17]. A general case of transforming arbitrary wavefronts using diffractive optical elements (DOEs) was also studied as a Monge-Kantorovich variational problem in a recent paper [18]. A proper ray mapping may be firstly computed by solving a minimization problem and then followed by finding an DOE to realize it. However, the numerical algorithms for this optimization problem were not detailed in their paper.

In our previous work [19], we simplify the problem by breaking up the problem into two separate steps too. Firstly, we compute the ray mapping based on Energy conservation using the variable separation method [7]. Then, we construct the two freeform optical surfaces simultaneously and point-by-point based on Snell’s law and the constancy of the optical path lengths (OPLs) according to this mapping. Such a simultaneous point-by-point calculation, which is pioneered by Benitz *et al*. for powerfully coupling two input wavefronts into two output ones [20], offers an effective way for controlling both the irradiance and wavefront. However, this method is only applicable to separable beam irradiances. In addition, it may generate small surface errors since only the slope of a point is considered to compute the next point on the first freeform surface.

We extend our previous work to tackle non-separable beam irradiances. The detailed design outline is presented in section 2. We firstly compute an appropriate ray mapping utilizing an adaptive mesh method based on the L^{2} Monge-Kantorovich problem. Then, we develop a modified simultaneous point-by-point procedure to construct the two freeform surfaces with smaller errors. In section 3, we verify this method by showing a challenging design example of projecting a complex irradiance distribution and a non-ideal wavefront. Finally, a brief summary is given in section 4.

## 2. Design method

The design problem is shown schematically in Fig. 1. Consider a laser beam is traveling along the positive *z* direction. The input beam at *z* = 0 is supposed to have an irradiance distribution *I _{in}*(

*x*,

*y*) and the wavefront close to it is denoted by position vectors

**S**= (

*x*,

*y*,

*z*). The output beam irradiance at

*z*=

*d*is prescribed by

*I*(

_{out}*x'*,

*y'*) and the wavefront close to it is denoted by

**T**= (

*x'*,

*y'*,

*z'*).

*n*,

_{1}*n*and

_{0}*n*are set as the the refractive indices of three mediums divided by the two freeform surfaces. The design goal is to specify

_{2}**P**= (

*x*,

_{1}*y*,

_{1}*z*) and

_{1}**Q**= (

*x*,

_{2}*y*,

_{2}*z*) on the two freeform surfaces in order to realize the above transformations. The two steps of the design are shown in 2.1 and 2.2, respectively.

_{2}#### 2.1 Ray mapping computation

The ray mapping from the input beam to the desired output beam is mainly governed by Energy conservation:

For separable input and output irradiance distributions, the ray mapping can be simply obtained as*x' = f*(

*x*) as

*y'*=

*g*(

*y*) [7]. In a general case for non-separable irradiance distributions, the ray mapping must take the form:

*x'*=

*f*(

*x, y*) and

*y'*=

*g*(

*x*,

*y*). By introducing this ray mapping into the right side of Eq. (1), we can obtain Eq. (2):wherein

*J*(

*f*,

*g*) is the Jacobian matrix of the ray mapping. Equation (2) alone can’t uniquely specify the ray mapping and need to be intertwined with the geometry of the optical elements for a rigorous solution [14, 18]. However, this adds great complexity to the design especially for handling double freeform optical surfaces.

Several methods based on the *L ^{2}* Monge-Kantorovich problem can help to specify an appropriate ray mapping (see e.g. [21–23].). For example, Haker’s method [22] was utilized by Bäuerle

*et al.*for approximating an optimum ray mapping in their irradiance tailoring approach [12]. Here, our alternative way to derive a proper ray mapping is the parabolic Monge-Ampère (PMA) adaptive mesh method presented by Sulman

*et al.*recently [23].

The optimal mapping of the *L ^{2}* Monge-Kantorovich problem can be characterized as the gradient of a convex potential

*u*(

*x*,

*y*) as shown in Eq. (3) [24, 25]:

Neumann boundary conditions specified by Eq. (4) are enforced at boundary grid points. It is obvious that the steady-state solution *u _{∞}* (t→∞) solves Eq. (4). Thus, the ray mapping can be determined by taking the spatial gradient of

*u*.

_{∞}The PMA adaptive mesh method provides a fairly simple way for solving the approximated optimal ray mapping. It is also convenient for us to solve the inverse ray mapping when the output irradiance distribution is much more complicated than the input one. Once the ray mapping is known, the input ray sequences can be defined by the input wavefront position vectors **S** and the input ray vectors **In,** and the output ray sequences can be defined by the output wavefront position vectors **T** and the output ray vectors **Out** [19].

#### 2.2 Double freeform-surface construction

Now that the input and output ray sequences have been derived, we need to construct the two freeform optical surfaces to achieve that transformation.

A simultaneous point-by-point construction procedure was shown in our previous work [19]. From two starting points, the first freeform surface is constructed firstly along the *y*-direction (or *x*-direction) and then along the *x*-direction (or *y*-direction) by intersecting the next input ray with the tangent plane of the current point while the second surface is obtained by equaling the optical path lengths (OPLs) between the input and output wavefronts. The construction of the first freeform surface is very analogous to Euler’s method for solving ODEs. The difference is that it is a two-dimensional integration of a partial differential equation (PDEs) using the language of extrinsic differential geometry [1, 2]. Although this method is effective, it may generate small surface errors and sufficiently small steps must be taken. The performance could be improved by introducing a second-order scheme. One possibility is using the curvature of the first freeform surface [19]. However, it is complicated since the surface curvature may not be computed a prior. Elmer proposed a modified method for generating the cross section curve of a trough reflector to minimize the surface normal discrepancies [26]. This is analogous to the modified Euler’s method which is also a second-order scheme. We extended this approach in the simultaneous point-by-point construction of the two desired freeform optical surfaces. Figure 2 shows this method schematically.

We need to generate two initial curves on the two freeform optical surfaces. Assume we have obtained two points **P**_{i,1} and **Q**_{i,1} of the two initial curves, wherein *i* = 1,2,3,…or n-1. We can describe the procedure by computing the next points **P**_{i + 1,1} and **Q**_{i + 1,1}. We can firstly calculate the normal vector **N**_{i,1} based on Snell’s law so that the ray emitted from **P**_{i,1} can be redirected to **Q**_{i,1}. Next, we compute **P***'*_{i + 1,1}, an approximate value of **P**_{i + 1,1}, by intersecting the input ray vector **In**_{i + 1,1} with the tangent plane of **P**_{i,1}. Once we have **P***'*_{i + 1,1}, we can get it’s corresponding point **Q***'*_{i + 1,1} on the second surface immediately based on the constancy of the OPLs between the input and output wavefronts. Next, we can acquire the normal vector **N***'*_{i + 1,1} at **P***'*_{i + 1,1} so that the ray emitted from **P***'*_{i + 1,1} can reach **Q***'*_{i + 1,1}. Next, we need to compute an intermediate point **P**_{i + 1/2,1} by intersecting **In**_{i + 1/2,1} with the tangent plane of **P**_{i,1}, wherein **In**_{i + 1/2,1} is emitted from the midpoint of **S**_{i,1} and **S**_{i + 1,1}. After that, we can obtain **P**_{i + 1,1} by intersecting **In**_{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} using the constant OPL condition. We can generate the entire two initial curves by repeat the process described above.

Similarly, we can derive all the points **P**_{i,2} and **Q**_{i,2} of the next two curves from the two initial curves, wherein *i* = 1,2,3,…n. We firstly calculate **P***'*_{i,2} by intersecting **In**_{i,2} with the tangent plane of **P**_{i,1}. Next, we can get **Q***'*_{i,2} based on the constancy of the OPLs. After that, We can acquire **N***'*_{i,2} based on Snell’s law in vector form. Then, we obtain **P**_{i,2} from the knowledge of **N**_{i,1}, **N***'*_{i,2} and an intermediate point **P**_{i,3/2}. We repeat this process to generate all the curves and interpolate them into the entire two freeform surfaces. Figure 3 illustrates the flow diagram of this point-by-point procedure.

## 3. Design example

To test the functionality of the proposed method, we design a Galilean type dual-lens beam shaper in a setting sketched in Fig. 4. The output irradiance distribution is prescribed as the letters “NIO” on a square background with an 256 × 256 grid (see Fig. 5(a)). The desired output wavefront is set as an astigmatism aberration with a peak-valley value (PV) of 80*λ*, wherein the beam wavelength *λ* = 1064nm (see Fig. 5(b)). The input laser beam is supposed to have a Gaussian irradiance distribution (beam waist: 5mm) and a spherical wavefront with a PV of 100λ.

We implement the PMA method to solve the inverse ray mapping *x* = *f'*(*x',y'*) and *y* = *g'*(*x',y'*) since the output irradiance distribution is much more complicated than the input one. The final adaptive mesh of the input irradiance is shown in Fig. 6(a). A uniform Cartesian mesh of the output irradiance is also presented in Fig. 6(b) for better visualizing the Energy transfer. No surprisingly, the ray mapping is more strongly deformed around the regions corresponding to the edges of the letters where the target irradiance is essentially discontinuous.

According to the ray mapping, we then compute the two required freeform surfaces using the modified simultaneous point-by-point procedure described in 2.2. In the calculation, the wavefronts immediately behind the input plane and before the output plane can be simply obtained by multiplying the input wavefront and the output wavefront by 1/1.5083, respectively. They are directly used in the simultaneous point-by-point procedure for integrating the two freeform optical surfaces. Figure 7 shows the two resulting freeform optical surfaces. They are fully continuous although the target irradiance is discontinuous. The first freeform surface is neither concave nor convex, which may add difficulties in its fabrication. Fortunately, we may obtain concave or convex freeform optical surfaces by crossing the ray mapping.

The angular spectrum method [27] is used to simulate the designed beam shaper with a FFT (Fast Fourier Transform) size of 2048 × 2048, wherein the input light field is interpolated into an 512 × 512 grid for a more accurate simulation. For this case, the two lenses of the beam shaper can be considered as pure phase elements in the simulation, and their phase distributions can be computed mainly based on their thickness functions [27]. The resulting output wavefront is derived by unwarping firstly the columns and then the rows of the wrapped phase [28]. The performances of the resulting output irradiance *I _{S}* is quantified by the relative root-mean-square deviation (RRMSD) from the target irradiance, as shown in Eq. (6):

*z'*is the

_{s}*z*-coordinate value of the simulated output wavefront.

To make a comparison, we also simulate a beam shaper with two freeform surfaces constructed by the previous simultaneous point-by-point procedure. Figure 8 shows the simulation results, and it is observed that the modified simultaneous point-by-point procedure brings improvements in both the output irradiance and the output wavefront. It seems that the wavefront errors shown in Figs. 8(b) and 8(d) are directly linked to the desired wavefront shown in Fig. 5(b). This indicates that the surface errors may be closely related to the shape of the output wavefront. Table 1 shows the performance parameters. They are mainly influenced by the ray mapping errors, the surface errors, and the diffraction effect, etc. Compared with the previous procedure, the modified one decreases the RRMSD of the prescribed irradiance distribution by about 18.9% and the RMS error of the prescribed wavefront by about 25.7%. The performance gap could be further enlarged when using less grid points. It is necessary to point out that the fabrication errors may destroy the improvements a lot.

## 4. Conclusion

In conclusion, we have demonstrated a general method to design two freeform optical surfaces for tailoring both the beam irradiance and wavefront. This method is beyond the limitation of our previous work [19] that the beam irradiance must be separable. Firstly, we determine an appropriate ray mapping using the PMA method [23] in the L^{2} Monge-Kantorovich framework. Then, we integrate the two freeform optical surfaces using a modified simultaneous point-by-point procedure to reduce the surface errors. Simulations of a challenging example yielded promising results in terms of the desired irradiance and wavefront distributions.

## Acknowledgments

The study was sponsored by the National Natural Science Foundation of China (Grant No. 61178055 and Grant No. 51021064) and the “973” Major State Basic Research Project of China (No. 2011CB706701).

## References and links

**1. **W. A. Parkyn, “Illumination lenses designed by extrinsic differential geometry,” Proc. SPIE **3482**, 389–396 (1998). [CrossRef]

**2. **W. A. Parkyn, and Jr., “Illuminating lens designed by extrinsic differential geometry,” Teledyne Lighting and Display Products, Inc., US Patent 5924788 (1999).

**3. **B. Parkyn and D. Pelka, “Free form lenses designed by a pseudo-rectangular lawnmower algorithm,” Proc. SPIE **6338**, 633808 (2006). [CrossRef]

**4. **W. A. Parkyn and D. G. Pelka, “Free-form lenses for rectangular illumination zones,” Anthony, Inc., US Patent 7674019 (2010).

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

**6. **V. Oliker, “Geometric and variational methods in optical design of reflecting surfaces with prescribed illuminance properties,” Proc. SPIE **5942**, 594207 (2005). [CrossRef]

**7. **L. Wang, K. Y. Qian, and Y. Luo, “Discontinuous free-form lens design for prescribed irradiance,” Appl. Opt. **46**(18), 3716–3723 (2007). [CrossRef] [PubMed]

**8. **Y. Ding, X. Liu, Z. R. Zheng, and P. F. Gu, “Freeform LED lens for uniform illumination,” Opt. Express **16**(17), 12958–12966 (2008). [CrossRef] [PubMed]

**9. **Y. Luo, Z. Feng, Y. Han, and H. Li, “Design of compact and smooth free-form optical system with uniform illuminance for LED source,” Opt. Express **18**(9), 9055–9063 (2010). [CrossRef] [PubMed]

**10. **F. R. Fournier, W. J. Cassarly, and J. P. Rolland, “Fast freeform reflector generation using source-target maps,” Opt. Express **18**(5), 5295–5304 (2010). [CrossRef] [PubMed]

**11. **D. Michaelis, P. Schreiber, and A. Bräuer, “Cartesian oval representation of freeform optics in illumination systems,” Opt. Lett. **36**(6), 918–920 (2011). [CrossRef] [PubMed]

**12. **A. Bäuerle, A. Bruneton, R. Wester, J. Stollenwerk, and P. Loosen, “Algorithm for irradiance tailoring using multiple freeform optical surfaces,” Opt. Express **20**(13), 14477–14485 (2012). [CrossRef] [PubMed]

**13. **R. Wu, L. Xu, P. Liu, Y. Zhang, Z. Zheng, H. Li, and X. Liu, “Freeform illumination design: a nonlinear boundary problem for the elliptic Monge-Ampére equation,” Opt. Lett. **38**(2), 229–231 (2013). [CrossRef] [PubMed]

**14. **H. Ries, “Laser beam shaping by double tailoring,” Proc. SPIE **5876**, 587607 (2005). [CrossRef]

**15. **T. Glimm and V. Oliker, “Optical design of two-reflector systems, the Monge-Kantorovich mass transfer problem and Fermat’s principle,” Indiana Univ. Math. J. **53**(5), 1255–1277 (2004). [CrossRef]

**16. **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]

**17. **J. Rubinstein and G. Wolansky, “Intensity control with a free-form lens,” J. Opt. Soc. Am. A **24**(2), 463–469 (2007). [CrossRef] [PubMed]

**18. **V. Oliker, J. Rubinstein, and G. Wolansky, “Ray mapping and illumination control,” J. Photonics Energy **3**(1), 035599 (2013). [CrossRef]

**19. **Z. Feng, L. Huang, M. Gong, and G. Jin, “Beam shaping system design using double freeform optical surfaces,” Opt. Express **21**(12), 14728–14735 (2013). [CrossRef] [PubMed]

**20. **P. Benítez, J. C. Miñano, J. Blen, R. Mohedano, J. Chaves, O. Dross, M. Hernández, and W. Falicoff, “Simultaneous multiple surface optical design method in three dimensions,” Opt. Eng. **43**(7), 1489–1502 (2004). [CrossRef]

**21. **J. D. Benamou and Y. Brenier, “A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem,” Numer. Math. **84**(3), 375–393 (2000). [CrossRef]

**22. **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]

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

**24. **R. T. Rockafellar, “Characterization of the subdifferentials of convex functions,” Pac. J. Math. **17**(3), 497–510 (1966). [CrossRef]

**25. **R. J. McCann, “Existence and uniqueness of monotone measure-preserving maps,” Duke Math. J. **80**(2), 309–323 (1995). [CrossRef]

**26. **W. B. Elmer, *The Optical Design of Reflectors*, 2nd ed. (Wiley, 1980), Chap.4.

**27. **J. W. Goodman, *Introduction to Fourier Optics*, 3rd ed. (Roberts and Company, 2005).

**28. **M. Gdeisat and F. Lilley, “Two-dimensional phase unwrapping problem,” http://www.ljmu.ac.uk/GERI/CEORG_Docs/Two_Dimensional_Phase_Unwrapping_Final.pdf.