A beam transformer converts a given incident beam into a prescribed output beam with certain irradiance and phase (or wavefront) distributions [1]. It has various applications, including lithography, material processing, laser or LED projector, optical communications, and light detection and ranging (lidar). Refractive, reflective, and diffractive optical elements can be used for different configurations of beam transformers [1]. Here, we focus on the commonly used refractive or reflective beam transformers, where ray optics are usually used in the design process. The design problem is mainly governed by three types of equations [2]: the energy conservation within a bundle of rays, the ray-tracing equations dominated by Snell’s law in vector form, and the Malus–Dupin theorem that describes the equal optical path lengths (OPLs) between the input and output wavefronts. Moreover, surface continuity should be considered for fabrication issues. Generally, the design problem is very difficult to handle because the output beam to be produced could have arbitrary amplitude and phase distributions without any a priori simplifying assumptions. Direct formulation of the design problem by merging the three types of equations into only one could be impossible, especially for the case that the output phase cannot be represented analytically.

Traditional beam transformers are commonly used for the cases in which the input and output beam wavefronts are kept planar and the irradiance distributions are rotationally symmetric. In the 1960s, Frieden [3] and Kreuzer [4] independently introduced a pair of plano-aspherical lens systems for converting a TEM00 Gaussian laser beam into a circular uniform beam while the output wavefront is kept planar. In their methods, the solutions for the aspheric surfaces can be formulated as an integral equation. Rhodes and Shealy [5] derived a set of ordinary differential equations for achieving the same goal. Many later designs can transform the Gaussian laser beam into collimated output beams with Fermi–Dirac, super-Gaussian, super-Lorentzian, and the flattened Lorentzian distributions (see, e.g., Ref. [6]). Doskolovich et al. [7] proposed an axisymmetric refractive optical element design for a more general case in which the input and output beam wavefronts are rotationally nonplanar.

Freeform beam transformer designs without rotational symmetries are mostly restricted to paraxial approximations. In a double-mirror system design, Nemoto et al. [8] first determined the coordinate relationships, i.e., the ray mapping, between the input Gaussian beam and the rectangular uniform output beam by applying variable separation to the energy conservation, and then obtained the first mirror surface by applying stationary phase method to the Fresnel diffraction integral and acquired the second mirror surface by compensating for the phase distortion. Shealy and Chao [9] designed a two-mirror system for transforming an elliptical Gaussian beam into a rectangular uniform output beam while keeping the beam collimated. After obtaining an analytical expression of the ray mapping that is independent in two orthogonal directions, they integrated the two reflective freeform surfaces explicitly from two first-order partial differential equations (PDEs). In the design of a plano-freeform lens pair, Feng et al. [10] calculated the variable-separable ray mapping numerically and constructed the double freeform surfaces point by point following the ray mapping. Such a method was extended for tackling nonseparable beam shaping by computing the ray mapping with the help of the L2 optimal transport (OT) theory [11,12]. However, Bösel et al. demonstrated that the L2 OT ray mapping is only accurate for paraxial or small-angle approximation in the design of double freeform surfaces for collimated beam shaping [13].

Most of the existing methods applicable for nonparaxial cases are devised for planar or spherical beam wavefronts. For these cases, the design problem can be formulated as a nonlinear PDE of the Monge–Ampère (MA) type [14,15]. The mathematical derivation starts with expressing the target plane coordinates as functions of the first freeform surface and its gradient, where the second freeform surface can be eliminated by imposing the OPL constancy condition. The resulting ray-tracing equations are then imported into the energy conservation between the incident and outgoing beam irradiance distributions to obtain the final MA equation, where the surface integrability condition is imposed to ensure a smooth surface. Such a direct determination has the advantage that the design problem can be described by only one equation but at the cost of an extremely complicated and tedious derivation process. Alternatively, collimated beam shaping problems for reflective or refractive designs can be described using OT theories that acquire an integrable ray mapping by finding a solution to a linear programming problem [16–19]. Rubinstein and Wolansky [20] 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. Oliker et al. [21] designed the required pair of plano-freeform lenses using the supporting quadric method (SQM) combined with OT. Mingazov et al. [22] proposed a new version of the SQM for collimated beam shaping by translating it into the gradient method of maximizing a concave function. Doskolovich et al. [23] showed that the computation of an integral ray mapping is reduced to an OT problem with a nonquadratic cost function that can be solved by a linear assignment algorithm. Yadav et al. [24] proposed a generalized least squares method for solving either a single lens with double freeform surfaces or two separate plano-freeform lenses in collimated beam shaping. Wei et al. extended [25] the symplectic flow-mapping scheme for irradiance tailoring [26] to the design of double freeform surfaces for transforming an elliptic Gaussian beam into a convergent beam with a complex irradiance distribution.

There are very few design methods available for complex input and output wavefronts. Feng et al. [27] provided a ray-mapping method in designing double freeform surfaces that can produce two prescribed irradiance distributions on two successive target planes. However, they still employ the restrictive L2 OT ray mapping, which is generally not accurate for a nonparaxial case. Bösel et al. [28] showed the ability to generate an astigmatic wavefront while forming a complex irradiance distribution. Instead of formulating the problem into one PDE, they directly solve the three nonlinear PDEs for the first freeform surface and the projected mapping coordinates, respectively, and then calculate the second freeform surface based on the OPL constancy. Such a strategy needs more computational power, and the ability to generate a more general wavefront without an analytical expression is not demonstrated.

We propose a new iterative wavefront tailoring (IWT)-based method to tackle the freeform lens design problem for producing a prescribed output beam with a complex irradiance distribution and phase profile. A previous IWT method is developed for designing a freeform lens that can generate a prescribed irradiance distribution on a planar or curved target through the intermediate construction of an outgoing wavefront immediately behind the exit freeform surface [29,30]. This method can simplify the formula derivation process and flexibly generate a variety of freeform optical structures with high accuracy. The major contribution of this work consists in how to deal with the difficulties arising from the additional generation of a complex output phase profile, which could realize optical field control using freeform optics in a simple, flexible, and accurate way. The method is described in detail in Section 2. We verify the method’s effectiveness in Section 3, where we design a double freeform lens for generating an output beam with complex irradiance and phase distributions that can produce a different irradiance distribution at a certain distance behind the first output plane. A comparison with the L2 OT ray mapping method is also included in Section 3.

We consider here the design of a freeform lens with two freeform surfaces performing a one-to-one mapping of the input optical field into a prescribed output optical field, as sketched in Fig. 1. The amplitude and phase distributions of the input beam are denoted by A(u,v) and ϕ(u,v) respectively, where (u,v)∈Ω⊂R2. The amplitude and phase distributions of the output beam are denoted by A′(ξ,η) and ϕ′(ξ,η), respectively, where (ξ,η)∈Γ⊂R2. An arbitrary light ray that goes through point S=(u,v,h) at the input plane strikes the first freeform surface at point P=(x1,y1,z1) and then strikes the second freeform surface at point Q=(x2,y2,z2). The output light ray intersects with the output plane at point T=(ξ,η,d). The input wavefront located at a given distance from the input plane is described by position vector W˜=(x˜,y˜,z˜), and the output wavefront located at a given distance from the output plane is described by position vector W˜′=(x˜′,y˜′,z˜′). The refractive indices of the medium at the left side of the first freeform surface and the right side of the second freeform surface are denoted by n1, and the refractive index of the freeform lens is denoted by n2.

In the irradiance control IWT method, an outgoing wavefront located behind the exit surface of the freeform lens to be designed is adopted to help solve a ray mapping. For the problem of optical field control, we employ an intermediate wavefront that is next to the entrance freeform surface to be designed. This intermediate wavefront is described by position vector W=(x,y,w), as illustrated in Fig. 1. In the following, we will first establish the parametrized intermediate wavefront equation, and then describe the IWT procedure specified for optical field control. For the sake of completeness, we also describe the double-freeform surface reconstruction from a ray mapping.

Since we consider a lossless system, and there is no crossing of rays from the input plane to the target plane, the irradiance distributions on the two planes should satisfy the energy conservation in differential form, I(u,v)=I′(ξ(u,v),η(u,v))|ξuηv−ξvηu|,(1)where I(u,v)=A(u,v)2, I′(ξ,η)=A′(ξ,η)2, and the subscripts here denote partial derivatives.

We now describe the light propagation properties inside the lens with the help of the intermediate wavefront W=(x,y,w). According to Fermat’s principle, the first partial derivatives of w with respect to x and y can be described as wx=−(x2−x)/(z2−w),wy=−(y2−y)/(z2−w).(2)

Both x and y can be thought as functions of the independent variables u and v, and according to chain’s rule, wu=wxxu+wyyu,wv=wxxv+wyyv.(3)

Combining Eq. (2) and Eq. (3), we have {wx=(yvwu−yuwv)/γ=−(x2−x)/(z2−w),wy=(xuwv−xvwu)/γ=−(y2−y)/(z2−w).(4)

Herein, γ=xuyv−xvyu. From Eq. (4), we can express (x2,y2) as {x2=x+(z2−w)(yuwv−yvwu)/γ,y2=y+(z2−w)(xvwu−xuwv)/γ.(5)

We denote the unit outgoing ray vectors as O=(O1,O2,O3). In the framework of geometrical optics, the phase is equivalent to the rays, and the ray directions are determined by the first partial derivatives of the phase (see, e.g., Refs. [31,32]). Therefore, the relationship between O and the output phase ϕ′(ξ,η) is O1=ϕξ′,O2=ϕη′,O3=1−ϕξ′2−ϕη′.(6)

(ξ,η) and (x2,y2) can be linked through a collinear relationship, $$\begin{gathered} {ξ=x2+(d−z2)O1/O3=x2+s,η=y2+(d−z2)O2/O3=y2+t. \\ \end{gathered}$$(7)

Herein, s means the displacement between ξ and x2, and t means the displacement between η and y2. Inserting Eq. (5) into Eq. (7), we obtain {ξ=x+s+(z2−w)(yuwv−yvwu)/γ,η=y+t+(z2−w)(xvwu−xuwv)/γ.(8)

Generally, s, t, and also z2 are dependent on ξ and η. For a degenerate case in which the output beam wavefront is planar, both s and t are zeros, and the difficulty mainly lies in handling z2. Such a problem is similar to the irradiance control problem on a curved target where the z-values are not constant [30]. Here, we focus on a more complicated case where the output phase is difficult to be expressed analytically. In such a case, s and t, which are directly relevant with the phase gradient, are difficult to be expressed explicitly with ξ and η.

Our solution for simplifying the complicated optical field control problem is to retain s and t on the right side of Eq. (8) and consider them as functions of the independent variables of u and v. Equation (8) and its differentiation with respect to u and v are inserted into Eq. (1) to eliminate the two variables ξ and η. The resulting equation is a second-order PDE of MA type, wuuwvv−wuv2+A1wuu+A2wuv+A3wvv+A4=0,(9)where the coefficients Ai(i=1,2,3,4) are functions depending on x,y,w,s,t,z2, and their first partial derivatives, and the second partial derivatives of x and y (see Appendix A). Rewrite Eq. (8) as a vector-valued function, (ξ,η)=m(u,v), and then a nonlinear boundary condition to Eq. (9) can be specified as m(∂Ω)→∂Γ,(10)where ∂Ω and ∂Γ are the boundaries of Ω and Γ, respectively. You may notice that Eq. (9) is unable to be solved unless x,y,s,t, and z2 are given in advance. However, x,y,s,t, and z2 are coupled with the lens data. To solve this contradiction, we introduce the following iterative procedure.

The IWT procedure for optical field control is illustrated in Fig. 2, and its detailed description is included as follows.

Step 0: We first define an input ray sequence described by the unit input ray vector I(u,v) and an input wavefront W˜(u,v) from the given input beam information. To start the iterative process, we need to provide an initial estimate of the ray mapping, (ξ,η)=(ξ(u,v),η(u,v)).

Step 1: From the ray mapping, we acquire the corresponding data of the first partial derivatives of the output phase, ϕξ′(ξ(u,v),η(u,v)) and ϕη′(ξ(u,v),η(u,v)). Then, we immediately determine an output ray sequence O(ξ(u,v),η(u,v)) according to Eq. (6). From the output ray sequence, we can reconstruct an output wavefront W˜′(ξ(u,v),η(u,v)) based on least squares, where we employ a concise relationship between the wavefront points and the outgoing ray vectors [27], {(W˜i+1,j′−W˜i,j′)·(Oi+1,j+Oi,j)=0,(W˜i,j+1′−W˜i,j′)·(Oi,j+1+Oi,j)=0.(11)

Equation (11) means that the chord linking the adjacent wavefront points is perpendicular to the average of the outgoing vectors through the two points. Equation (11) can be written as a linear system of equations of the distances between T and W˜′, which can be solved with least squares based on Hermann’s method [33] and the MATLAB backslash operator [34].

Step 2: Since we have determined the input ray sequence and wavefront, and the output ray sequence and wavefront, we can compute the double freeform surfaces P(u,v)=(x1(u,v),y1(u,v),z1(u,v)) and Q(u,v)=(x2(u,v),y2(u,v),z2(u,v)) to realize the desired transformation using least squares based on Snell’s law and the condition of the OPL constancy [27]. More details are included in Subsection 2.C.

Step 3: After obtaining the data of the double freeform surfaces, we can determine s and t. We can also obtain an intermediate ray sequence R=Unit(Q−P)=(Q−P)/|Q−P|, from which we reconstruct an intermediate wavefront W=(x,y,w˜) using least squares. Herein, w˜ is used to differentiate from w in Eq. (9).

Step 4: We insert the data of x,y,s,t, and z2 into Eq. (9) and Eq. (10). Taking w˜ as the initial value, we solve the MA equation numerically. After obtaining a solution w, we can acquire a new ray mapping, (ξ,η)=(ξ(u,v),η(u,v)), based on Eq. (8).

Step 5: The new ray mapping is inserted into Step 1 to repeat the loop. The stop criterion can be chosen as the difference between the current ray mapping and the previous one or a certain number of iterations.

Step 6: After obtaining the satisfied ray mapping, we need to specify its corresponding output ray sequence and wavefront as in Step 1 and reconstruct the final double freeform surfaces as in Step 2.

Such an iterative procedure can make the freeform lens design for complicated optical field control easier to implement. Although a sequence of MA equations needs to be solved, a multiscale strategy [29,30] for the irradiance control IWT methods can be adapted to speed up the computation here.

Here, we recall the reconstruction of double freeform surfaces following the one-to-one mapping between an input ray sequence and an output ray sequence using a least squares-based iterative procedure [27], which is effective and fast. The basic steps are as follows.

Step 0: We first give an initial approximation of the entrance freeform surface, which can be simply set as a plane.

Step 1: After that, we can obtain a corresponding exit freeform surface based on the constant OPL condition, n1[W˜,P]+n2[P,Q]+n1[Q,W˜′]=Const,(12)where [A,B] refers to the distance between A and B.

Step 2: The required unit ray sequence R inside the lens can be acquired as R=Unit(Q−P). After that, the required normal field of the first freeform surface can be acquired based on Snell’s law in vector form, N=Unit(n2R−n1I).

Step 3: Notice that the normal field N is not necessarily integrable and there generally does not exist an exact freeform surface that is perpendicular to N everywhere. Therefore, we also employ a least squares method to compute an approximate freeform surface P, and the basic relationship between P and N is [27] {(Pi+1,j−Pi,j)·(Ni+1,j+Ni,j)=0,(Pi,j+1−Pi,j)·(Ni,j+1+Ni,j)=0.(13)

Equation (13) can be written as a linear system of equations of the distances between W˜ and P, which can also be solved using Hermann’s method and the MATLAB backslash operator.

Step 4: The updated entrance freeform surface P is inserted into Step 1 to obtain a new exit freeform surface Q. The above process is repeated until a required stop criterion is satisfied.

As a demonstrative example, we design a freeform lens with the refractive index of 1.5 for transforming a divergent elliptical Gaussian beam into an specified output beam with a complex irradiance distribution and phase profile (see Fig. 3). The input beam is emitted from a point light source located at the original point, and the 1/e2 full width divergence angles for x and y directions are 12° and 36°, respectively. Figure 3(a) shows the input beam irradiance I(u,v) on the z=20mm plane, where (u,v)∈Ω={(u,v)||u|≤2.1021mm,|v|≤6.4984mm}. The input beam phase distribution ϕ(u,v), which corresponds to a spherical wavefront, is shown in Fig. 3(b). Figure 3(c) shows the desired output beam irradiance distribution I′(ξ,η) on the plane of z=45mm, where (ξ,η)∈Γ={(ξ,η)||ξ|≤12mm,|η|≤12mm}. Figure 3(d) shows the desired output phase distribution ϕ′(ξ,η), which is rather irregular and can be difficult to describe by an analytical formula. This phase distribution is aimed for converting I′(ξ,η) into a different irradiance distribution I′′(ξ′,η′) on the plane of z=75mm [see Fig. 3(e)], where (ξ′,η′)∈Γ′={(ξ′,η′)||ξ′|≤16mm,|η′|≤16mm}. ϕ′(ξ,η) is obtained by numerically solving a phase-retrieval problem from the two output irradiance distributions [35], $$\begin{gathered} ϕξξ′ϕηη′−ϕξη′2+χ(1−ϕη′2)ϕξξ′+2χϕξ′ϕη′ϕξη′+χ(1−ϕξ′2)ϕηη′ \\ +χ4d′2{1−I′(ξ,η)/I′′[(ξ,η)+∇ϕ′/χ]}=0, \end{gathered}$$(14)where d′ denotes the distance between the two output planes, χ=1−|∇ϕ′|2/d′, and the coordinates (ξ′,η′) for describing I′′(ξ′,η′) have been replaced with (ξ,η)+∇ϕ′/χ.

The computation size is desired as 256×256. To speed up the computation, we employ a multiscale strategy that is similar to that in Refs. [29,30]. The initial computation size is set to 32×32. In this step, the initial ray mapping is obtained from solving a standard MA equation corresponding to the L2 OT (see e.g., Ref. [27]), ψuuψvv−ψuv2=I(u,v)/I′(∇ψ),(15)where ψ is a potential whose gradient ∇ψ=(ξ,η). After implementing the computations on the grid of 32×32 points for three iterations, the final 32×32 points ray mapping is linearly interpolated into the one of 64×64 points, which is used to start the iterative computations on the grid of 64 × 64 points. Such a process is repeated until implementing the computations on the grid of 256×256 points for three iterations. The final ray mapping specifies the final output ray sequence Oi,j and wavefront W˜i,j′, both of which are used for obtaining the final entrance and exit freeform surfaces according to Section 2.C. We solve the standard MA equation shown in Eq. (15) and the generalized MA equation shown in Eq. (9) [and also Eq. (14)] following the numerical technique outlined in Ref. [29]. A Newton–Krylov solver, Kelley’s nsoli.m is adopted to solve the nonlinear system of equations obtained from discretizing the MA equations [36], where maxit (the maxmium number of nonlinear iterations) is set as 200, maxitl (the upper bound for linear iteration per nonlinear iteration) is set as 20, and etamax (the maximum error tolerance for residual in inner iteration) is set as 0.1. Since I′(ξ,η) can be difficult to represent analytically, we employ linear interpolation to acquire its values corresponding to the scattered grid points of (ξ,η) in each iteration of the Newton–Krylov solver. We also use linear interpolation to obtain ϕξ′(ξ,η) and ϕη′(ξ,η) when obtaining a ray mapping in each IWT iteration. The computations are implemented in MATLAB 2020a on an Intel Core i7-6700 k CPU at 4.0 GHz with 48 GB RAM, running Windows 10 (64 bit), which takes around 2.43 s, 11.68 s, 40.11 s, and 93.93 s for the grid points of 32×32, 64×64, 128×128, and 256×256, respectively. The solution to the MA equations accounts for more than 90% of the total running time.

Let ϵ denote the mean absolute deviation of the angles between the computed normals to the entrance freeform surface, i.e., Pu×Pv, and the required normals for reconstructing the entrance freeform surface in each IWT iteration. Figure 4 illustrates the effectiveness of the proposed multiscale IWT algorithm in the reduction of the ϵ value. The ϵ value for the final entrance freeform surface is 0.0064°. Figures 5(a) and 5(b) show the final entrance and exit freeform surfaces, respectively, where you can observe structures that can cause difficulties in analytical expression.

The scattered data points of the double freeform surfaces are inserted into Rhinoceros to create a solid lens model. We implement Monte Carlo ray tracing with 2×107 rays in LightTools 9.0 to demonstrate the performance of the designed freeform lens. Figure 6(a) visualizes the ray-tracing results of the designed freeform lens in LightTools. Figures 6(b) and 6(c) show the simulated irradiance distributions on the first and second output planes, respectively. It can be seen that the simulated irradiance distribution on the first output plane agrees well with the desired output irradiance distribution shown in Fig. 3(c). The simulated irradiance distribution on the second output plane is also in line with the desired one shown in Fig. 3(e), implying that the output phase on the first output plane is achieved well.

We provide an L2 OT ray-mapping design for comparison. The L2 OT ray mapping is also calculated in a multiscale way: the computed potential ψ by solving Eq. (15) for a coarser grid is linearly interpolated into an initial value for a finer grid. The double-surface reconstruction following the L2 OT ray mapping is the same as that described in Section 2.C. Since the design geometry is far beyond the paraxial and small-angle approximation, the ϵ value for the L2 OT ray-mapping design is as high as 2.3197°, resulting in severe deformations in the simulation results, as shown in Figs. 6(d), 6(e), and 6(f).

We employ the relative root-mean-square deviation (RRMSD) to quantify the deviation of the simulated irradiance distribution from the prescribed one: RRMSD=∥Ipres−Isimu∥F/∥It∥F, where Ipres and Isimu denote the matrices of the prescribed and simulated irradiance distributions, respectively, within the given domain, and ∥·∥F refers to the Frobenius norm [37]. For the L2 OT ray-mapping method, the RRMSD values for the two simulated irradiance distributions can be higher than 0.329. The RRMSD values can be greatly reduced to be lower than 0.045 when we employ the new proposed method.

Figure 7(a) provides the final ray mapping for the new proposed method. Figure 7(b) shows the L2 OT ray mapping. It seems that the two ray mappings demonstrate similar deformations: those regions corresponding to the edges of the numerals are more strongly deformed. However, sufficient differences exist between the two ray mappings, as visualized in Fig. 7(c). The average distance between the two ray mapping is 0.2775 mm. As can be seen from Fig. 7(c), the ray-mapping differences around the central region have relatively small amplitudes, demonstrating that L2 OT ray mapping is more accurate for the paraxial region. However, it does not seem that the central regions on the two simulated irradiance distributions have better performances than the surrounding regions. The holes in the two simulated irradiance distributions of the L2 OT ray-mapping design are caused by a small “bump” around the center of the entrance freeform surface. This small bump is a result of implementing least squares surface reconstruction with a preset value of the center surface point when the normal error is large. We can observe whirls in Fig. 7(c), which may reflect the fact that the L2 OT ray mapping cannot lead to curl-free surface integration. Note that larger errors will occur when we reconstruct the double freeform surfaces according to the L2 OT ray mapping using a point-by-point integration strategy, as described in Ref. [11].

The discrepancies between the new proposed method and the L2 OT ray-mapping method can be narrowed for those designs which satisfy or do not deviate much from paraxial or small-angle approximation. Figure 8 provides a comparison of two freeform lenses designed with the two methods, respectively, for a case that the input beam becomes circular (the 1/e2 full width divergence angles for both x and y directions are 36°). The RRMSD values for the simulated irradiance distributions of the two freeform lenses become very close. The reason is that the L2 OT ray mapping becomes much closer to the ray mapping computed by the new proposed method. The average distance between the two ray mappings is 0.0275 mm. However, the freeform lens designed with the new proposed method still performs better enough in terms of surface accuracy. The ϵ value for the freeform lens designed with the new proposed method is 0.0071°, while the ϵ value for the freeform lens designed with the L2 OT ray-mapping method is 0.2119°. The normal error difference between the two lenses can result in larger irradiance discrepancies on receiving planes far away from the lenses.

We have proposed a new IWT-based method in designing double freeform optical surfaces for transforming a given optical field into a prescribed one. This method simplifies the design by gradually tailoring an intermediate wavefront next to the entrance freeform surface and incorporating the data of the output phase gradient and the z coordinate of the second freeform surface associated with the ray mapping into the IWT procedure. The solution to the wavefront equation leads to a ray mapping between the input beam and the output beam. Approximate first and second freeform surfaces are reconstructed following the ray mapping via an iterative reconstruction procedure based on ray-tracing equations and the constancy of OPLs. The new reconstructed double freeform surfaces can be used to obtain an updated wavefront equation that brings a new ray mapping. Such an iterative process can be implemented in a multiscale way to reduce calculation complexities. Results show that we can obtain a double-freeform lens generating two different irradiance distributions on two successive target planes with much better performances than those of the L2 OT ray mapping method for a case that deviates much from the paraxial or small-angle approximation.

Compared with the previous IWT methods for irradiance control, where a single freeform surface could be enough [29,30], the proposed method can realize additional generation of a prescribed output phase distribution requiring at least double freeform surfaces. This method is also applicable for designing double freeform reflective surfaces or a double plano-freeform lenses system. Such a special iteration strategy may also work with other methods, including the direct determination.

Currently, this method is restricted to those cases in which there are no zero values in the output beam irradiance distributions and there are no extreme local curvatures (or singularities) in the lens surfaces and the input and output wavefronts. We also need to mention that the starting geometry of the lens should be chosen to avoid unphysical solutions.