1. INTRODUCTION

In recent years, the manufacturing of freeform surfaces has become increasingly feasible. These freeform surfaces offer an elegant way of simultaneous irradiance and phase control. Therefore, the development of numerical algorithms for the calculation of continuous freeform surfaces for control of the irradiance and/or the phase of a beam is of great interest.

In this work, the problem of designing continuous double freeform surfaces in a geometrical optics approximation is considered, in which two collimated beams of arbitrary irradiance are mapped onto each other. The simultaneous control of the phase thereby requires the second freeform surface. Hence, in contrast to the single freeform surface design for irradiance control, additional difficulties arise, since the surfaces can in general not be calculated independently. Several methods for phase and irradiance control with double freeform surfaces have been proposed in the literature.

One of the first design methods for the control of two sets of wavefronts by coupled freeform surfaces is the simultaneous multiple surface (SMS) method, which was developed by Benitez and Miñano [1]. The surfaces are thereby constructed from generalized Cartesian ovals and by applying constant optical path length (OPL) conditions [2]. The design method can be utilized for numerous applications in imaging and nonimaging optics [2–4].

In Ref. [5], the author reports to have solved the problem of laser beam shaping with double freeform surfaces by solving a Monge–Ampère-type equation. Details about the algorithm are not given.

Zhang et al. [6] and Chang et al. [7] solve the design problem by describing it in the form of a Monge–Ampère-type partial differential equation (PDE), discretizing the equation by finite differences and then solving the resulting nonlinear equation system by the Newton method. The design method can be applied to a variety of wavefront shapes [7]. Limitations and the performance of the method are unfortunately not discussed.

An alternative approach to construct freeform surfaces for irradiance and phase control is from ray-mapping methods [8–19]. These methods are based on the separation of the design process into two separate steps: the calculation of an integrable ray mapping between the source and the target irradiance, and the subsequent construction of the continuous freeform surfaces from the mapping. The integrability thereby ensures the continuity of the freeform surfaces and the mapping of the input irradiance onto the ouptut irradiance. Since the integrability depends on the physical properties of the optical system, it is in general a nontrivial task to find such a mapping.

As shown in several publications, there is a strong relation between the inverse problem of nonimaging optics and optimal mass transport (OMT) [20–25]. The cost function, which has to be applied to a certain optical configuration, is thereby problem-specific. For example, the mapping of two collimated beams with arbitrary irradiance onto each other with double freeform mirrors is described by a quadratic cost function [21] and can be solved by corresponding numerical schemes [19]. The same problem statement with double freeform lenses, which is considered in this work, is described by a different cost function, which depends on the OPL between the freeform surfaces as it was shown by Rubinstein and Wolansky [23].

The investigations presented here are inspired by several publications [9–14], in which the authors applied the quadratic OMT cost function to calculate a ray mapping to deal with the lens design problem. With this ray mapping, designs have been demonstrated of both single freeform surfaces for irradiance control for collimated input beams and point sources [9,10,14], and that of double freeform surfaces for irradiance and phase control [11–13]. As demonstrated for illumination control with single freeform surfaces in Ref. [17] and for collimated beam shaping with double freeform surfaces in Ref. [18], the design problems are thereby restricted to a paraxial approximation. Design methods for irradiance and phase control with continuous double freeform lenses using the quadratic cost function ray mapping like [11,12,18] are therefore inherently restricted to the paraxial regime. This means that either large enough distances between the freeform surfaces have to be considered or the surface continuity has to be given up, which leads to more difficult manufacturing processes and possibly to stronger diffraction effects. The main purpose of this work is therefore to present a detailed and efficient numerical algorithm to overcome these restrictions to the paraxial regime. Additionally, it is argued theoretically that the convergence of the presented algorithm can be ensured, which is of key importance for freeform design methods dealing with highly nonlinear PDEs.

Therefore, we investigate the design by ray-mapping methods of double freeform surfaces that map two collimated beams with arbitrary irradiance onto each other beyond the paraxial approximation. To overcome the restriction to the paraxial regime, which is necessary for the construction of compact systems, the design problem will be modeled by two coupled PDEs. This involves on one hand the Jacobian equation, expressing the local energy conservation, and on the other hand a ray-mapping PDE, enforcing the surface continuity and the constant OPL. The PDEs will then be solved by a root-finding scheme with the OMT mapping from the quadratic cost function as the initial iterate, leading to a construction approach for the freeform surfaces.

To do so, the work is structured as follows. In Section 2, by using the law of refraction, the constant OPL condition and a surface continuity condition, a PDE for an integrable ray mapping, is derived. Together with the Jacobian equation it builds a system of PDEs for the determination of the mapping components. It is argued that the PDE system is fulfilled within the paraxial approximation by the quadratic cost function OMT map. In Section 3, a method for solving the PDE system for general distances between the freeform surfaces is presented. It is based on discretizing the PDEs with finite differences and solving the resulting system of nonlinear equations by a standard root-finding scheme with the quadratic cost function OMT map as the initial iterate. A summary of the design algorithm and a detailed discussion of the implemenation is presented in Section 4, followed by the application of the presented method to the design of a single-lens and a two-lens system in Section 5. Finally, in Section 6, a short discussion of the results is presented.

2. FREEFORM DESIGN IN PARAXIAL APPROXIMATION

A. Energy Conservation and Cost Functions

In Ref. [17], a design method was presented for the construction of a single freeform surface for a collimated input beam with irradiance $I_S(x,y)$ and an arbitrary illumination pattern $I_T(x,y)$ on a target plane. It was shown that in the paraxial approximation, the design process can be decoupled into two separate steps. In the first step, a ray mapping $\mathbf{u}(x,y)=(u_x(x,y),u_y(x,y))$ is calculated from the theory of OMT, and in the second step the freeform surface is constructed from the mapping.

There are several basic physical principles that a ray mapping needs to fulfill. First, to map the source irradiance $I_S(x,y)$ onto the target irradiance $I_T(x,y)$, the ray mapping should be energy conserving. The local energy conservation is expressed through the Jacobian equation

As shown by Rubinstein and Wolansky, the cost function for collimated beam shaping with double freeform lenses takes a different form than Eq. (2) [23]. The authors propose to minimize the corresponding cost function by a steepest-descent algorithm to get the ray mapping [23], but unfortunately, a numerically stable implementation is a nontrivial problem.

Due to its applicability in the paraxial approximation (see below) and the availability of numerous published stable numerical schemes for its calculation, the quadratic cost function OMT mapping serves as an initial iterate for the root-finding scheme presented below. It will therefore build the basis of the design approach presented in Sections 3 and 4.

B. Ray-Mapping Condition

We follow the approach from Refs. [17,18] by expressing the basic geometry according to Fig. 1 in terms of the collimated input and output vector fields ${\mathbf{s}}_1$ and ${\mathbf{s}}_3$, the refracted vector field ${\mathbf{s}}_2$ and the ray-mapping vector field ${\mathbf{s}}_4$

 

Fig. 1. (a) Irradiance distributions $I_S(x,y)$ and $I_T(x,y)$ with the boundaries $\partial {\mathrm{\Omega }}_S$ and $\partial {\mathrm{\Omega }}_T$ are given on the planes $z=0$ and $z=z_T$, respectively. In the first step an integrable ray mapping $\mathbf{u}(x,y)=(u_x(x,y),u_y(x,y))$ is calculated between the distributions, which defines the vector field ${\mathbf{s}}_4$ between source and target plane. (b) In the second step the freeform surfaces $z_{\mathrm{I}}(x,y)$ and $z_{\mathrm{II}}(x,y)$ are calculated from the law of refraction and the constant OPL condition. The collimated input and output beams are represented by the vector fields ${\mathbf{s}}_1$ and ${\mathbf{s}}_3$, and the refracted beam by the vector field ${\mathbf{s}}_2$ [18].

Download Full Size | PDF

Since the goal is to calculate at least continuous freeform surfaces, we have to apply the surface continuity condition

The left-hand side (LHS) of Eq. (8) represents the dot product of the projected gradient of the first surface $\nabla z_{\mathrm{I}}(x,y)=({\partial }_x{z}_{\mathrm{I}}(x,y),{\partial }_y{z}_{\mathrm{I}}(x,y))$ and the direction perpendicular to the ray deflection $(\mathbf{u}(x,y)-\mathbf{Id})$. Therefore, a nonvanishing RHS of Eq. (8) contradicts the law of refraction, which states that the incoming, the refracted, and the normal vector have to lie in the same plane. This can be seen directly by using the relation

A further relation can be derived by applying the chain rule to the equation ${\nabla }_{\mathbf{u}}{z}_{\mathrm{II}}(u_x,u_y)=({\partial }_{u_x}{z}_{\mathrm{II}}(u_x,u_y),{\partial }_{u_y}{z}_{\mathrm{II}}(u_x,u_y))=\frac{{\mathbf{n}}_{\mathrm{II}}(u_x,u_y)}{{({\mathbf{n}}_{\mathrm{II}})}_z}$ [27]. This provides us with the gradient $\nabla z_{\mathrm{II}}(x,y)$, which is used to rewrite the second term of the RHS of Eq. (8).

Hence from the continuity condition [Eq. (5)] and the law of refraction [Eq. (6)] follow the system of equations

These Eqs. (10) and (11), together with Eq. (1), build a system of PDEs for the unknown mapping $\mathbf{u}(x,y)$ and the surfaces $z_{\mathrm{I}}(x,y)$ and $z_{\mathrm{II}}(x,y)$.

To decouple the design process into separate steps as described in the beginning of this section, one needs to find a ray mapping fulfilling the condition [Eq. (11)] that is nontrivial without any a priori knowledge about the freeform surfaces. For single and double freeforms, this can be done by considering the small-angle approximation $(z_{\mathrm{II}}(u_x,u_y)-z_{\mathrm{I}}(x,y))\gg |\mathbf{u}(x,y)-\mathbf{Id}|$ leading to a vanishing first and second term in Eq. (11), since in this limiting case we get [using Eq. (4)] ${\mathbf{n}}_{\mathrm{I}}·{\mathbf{s}}_2\sim (z_{\mathrm{II}}(u_x,u_y)-z_{\mathrm{I}}(x,y))$, ${({\mathbf{n}}_{\mathrm{II}})}_z·|{\mathbf{s}}_2|\sim (z_{\mathrm{II}}(u_x,u_y)-z_{\mathrm{I}}(x,y))$, and the numerators are only dependent on the map $\mathbf{u}(x,y)$. Additionally, using the mapping from the quadratic cost function defined by Eq. (3), which gives $\nabla \mathbf{v}=0$, the condition [Eq. (11)] is fulfilled. Hence, the surfaces can be calculated by using Eq. (10) with appropriate boundary conditions. The boundary conditions can be derived by considering the law of refraction [Eq. (9)] on the boundaries of $I_S(x,y)$ and $I_T(x,y)$ [17]. As discussed in Ref. [17], this leads to a path-independent integration of Eq. (9) to calculate the surface $z_{\mathrm{I}}(x,y)$.

An alternative way to derive the validity of the quadratic cost function in the paraxial approximation results is by utilizing an expansion of the Rubinstein–Wolansky cost function [23] for small angles.

The double freeform design process can further be simplified by also using the constant OPL condition

Since ${\mathbf{n}}_{\mathrm{I}}·{\mathbf{s}}_2$ and ${({\mathbf{n}}_{\mathrm{II}})}_z·|{\mathbf{s}}_2|$ depend on $z_{\mathrm{II}}(u_x,u_y)-z_{\mathrm{I}}(x,y)$, Eq. (11) can be written as a PDE for the components of mapping function $\mathbf{u}(x,y)$ only. Therefore, Eqs. (1) and (11) build a system of PDEs for the functions $u_x(x,y)$ and $u_y(x,y)$. Both equations build the basis of the design process for double freeform surfaces described in Sections 3 and 4.

Before we give any details, we want to discuss the condition [Eq. (11)] briefly for freeform mirrors.

1. Freeform Mirrors

For mirrors, the refractive indices in Eq. (6) have to be replaced by $n_1\equiv n_2\equiv -1$, and we get ${\mathbf{n}}_{\mathrm{I}}{\mathbf{s}}_2=-({\mathbf{n}}_{\mathrm{II}})|{\mathbf{s}}_2|$. Therefore, Eq. (11) reduces to

Hence, in contrast to the single-lens, single-mirror and double-freeform lens systems, the design problem can be solved by the quadratic cost function without any additional assumptions like the paraxial approximation.

3. FREEFORM LENS DESIGN BEYOND PARAXIAL APPROXIMATION

As mentioned in the previous section, Eqs. (1) and (11) are the basis of the design approach presented in the following. Since Eq. (11) is exactly fulfilled by the mapping with the quadratic cost functions for an infinite distance between the freeform surfaces or an infinite OPL, respectively, Eq. (11) represents a correction to Eq. (3) for finite OPLs. Hence, for finite distances between $z_{\mathrm{I}}(x,y)$ and $z_{\mathrm{II}}(x,y)$, we are searching for corrections $\mathrm{\Delta }\mathbf{u}(x,y)=(\mathrm{\Delta }{u}_x(x,y),\mathrm{\Delta }{u}_y(x,y))$ with $\mathrm{\Delta }\mathbf{u}(x,y)\stackrel{\mathrm{OPL}\to \infty }{\to }0$ to the ray mapping defined by the quadratic cost function, which we will denote by ${\mathbf{u}}^{\infty }(x,y)$ in the following. Hence, after writing Eq. (11) in terms of Eq. (3) plus a perturbation term, we want to solve the system of equations

The scalability of the mapping correction $\mathrm{\Delta }\mathbf{u}(x,y)$ with the parameter ${\mathrm{OPL}}_{\text{red}}$ thereby suggests solving Eqs. (15) and (16) by a root-finding method, since ${\mathrm{OPL}}_{\text{red}}$ can be reduced step by step with the solution $\mathrm{\Delta }\mathbf{u}(x,y)$ from the previous step as an initial iterate. This ensures the convergence of the root finding and can be done until the desired design goal or ${\mathrm{OPL}}_{\text{red}}$, respectively, is reached.

Additionally, we want to apply boundary conditions to Eqs. (15) and (16). To do so, we use standard boundary conditions for OMT problems by demanding that the boundary of the support of $I_S(x,y)$ is mapped onto the boundary of the support of $I_T(x,y)$. In the case of mapping two unit squares onto each other, such as in the example in Section 5, we therefore apply

We discretize Eqs. (15) and (16) using the standard central finite differences for the derivatives of $\mathrm{\Delta }{u}_x(x,y)\to {(\mathrm{\Delta }{u}_x)}_{i;j}$ and $\mathrm{\Delta }{u}_y(x,y)\to {(\mathrm{\Delta }{u}_y)}_{i;j}$ with $i=1,\dots ,N$; $j=1,\dots ,N$ at the inner points $i=2,\dots ,N-1$; $j=2,\dots ,N-1$ and second-order finite differences at the boundary. For the derivatives of, e.g., $\mathrm{\Delta }{u}_x(x,y)$, we get

4. DESIGN ALGORITHM

The design algorithm for the construction of double freeform surfaces presented in the previous section is summarized in the following. Since some of the steps offer freedom to choose between different applicable numerical techniques, we add some remarks below which are important for the examples in Section 5.

There are numerous publications presenting numerical methods for the calculation of OMT maps for quadratic cost functions ${\mathbf{u}}^{\infty }(x,y)$. In the example section below, we choose the numerical scheme developed by Sulman et al. [28]. It provides an efficient calculation of ${\mathbf{u}}^{\infty }(x,y)$, but has some drawbacks, such as the restriction to square supports of the irradiance distributions $I_S(x,y)$ and $I_T(x,y)$. In our implementation of Sulman’s algorithm, we recognized instabilities if the distributions show large irradiance gradients. For such distributions in the example in Section 5 we therefore choose a background irradiance $I_T(x,y)>0$ to ensure the convergence. An OMT method that is less restrictive can be found, for example, in Ref. [29].

For solving the system of nonlinear equations in step 3 we have used the fsolve() function from MATLAB’s 2015b optimization toolbox with the trust-region-reflective solver. Providing fsolve() with the structure of the Jacobian matrix of the object function allows an efficient calculation of the solution of the nonlinear equation system even for a large number of variables. The scalability with the parameter ${\mathrm{OPL}}_{\text{red}}$ of the distance of the initial iterate to the solution of Eqs. (15) and (16) suggests that the root finding could be accelerated by using, e.g., the Newton algorithm. In practice, we never experienced the necessity of iteratively reducing ${\mathrm{OPL}}_{\text{red}}$ to ensure the convergence of the algorithm. It was always possible to directly set ${\mathrm{OPL}}_{\text{red}}$ to the final design goal.

Solving the equation system in the form of Eqs. (15) and (16), we recognized oscillations in the solution $\mathrm{\Delta }\mathbf{u}(x,y)$. These seem to arise due to the irradiance root finding by the Jacobian Eq. (15), which in contrast to Eq. (16), is already fulfilled to a high degree by the initial iterate $\mathrm{\Delta }\mathbf{u}(x,y)=0$. Hence, we replace the $\det (\nabla {\mathbf{u}}^{\infty }(x,y))I_T(\mathbf{u}(x,y))$ term in Eq. (15) by

To calculate the first freeform, we use Eq. (9) with Eq. (13), which can be integrated along an arbitrary path to give the surface sag. The basic assumption for this path-independent integration is the satisfaction of Eq. (16), as mentioned in Section 2. Since Eq. (16) cannot be perfectly fulfilled numerically, the path integration leads to an accumulation of errors and therefore to deviations from the predefined plane wavefront.

For the second surface, one could use Eq. (13) directly, which gives this surface on a scattered grid. Since we want to test the design algorithm with ray-tracing software and due to advantages for possible freeform manufacturing processes, we calculate the second surface on a regular grid. This can either be done by scattered data interpolation, which might lead to inaccuracies in the irradiance pattern, or by applying the ray-tracing formula

5. EXAMPLES

In the following, the presented design algorithm is applied to two example target distributions. The first example consists of redistributing a Gaussian input beam with a waist of $1/\sqrt{2}$a.u. onto the letters “IAP,” and in the second example, the same Gaussian will be mapped onto the test image “house” (see Fig. 2). In “IAP,” the difficulties arise mainly due to the steep gradients, whereas “house” shows numerous gray levels. For the resolution of both images, $250\mathrm{pixels}\times 250\mathrm{pixels}$ are chosen.

 

Fig. 2. (a) Input irradiance $I_S(x,y)$; (b) output irradiance $I_T(x,y)$ of the example “IAP;” (c) output irradiance $I_T(x,y)$ of the example “house.” The irradiances are normalized to ensure energy conservation.

Download Full Size | PDF

Since input and on unit squares are of the same size, the freeform surfaces have the shape of squares with the side length 1a.u. For the example “IAP,” we choose as the physical parameters of the system, the refractive indices $n_1=1.5$ of the lenses, $n_2=1$ of the surrounding medium and the desired reduced OPL ${\mathrm{OPL}}_{\text{red}}=-0.2$a.u. are used. Therefore, a two-lens system, such as in Fig. 1(b) is considered. A distance between both freeform surfaces of approximately 0.4a.u. can be estimated from Eq. (13), since, due to symmetry reasons, the corner points of the irradiances are mapped onto each other. The reference ray for calculating the OPD intersects $z_{\mathrm{I}}(-0.5,-0.5)$ and $z_{\mathrm{II}}(-0.5,-0.5)$.

For the example “house,” we calculate a single-lens system defined by the parameters $n_1=1$, $n_2=1.5$, and ${\mathrm{OPL}}_{\text{red}}=0.2$a.u., leading to an approximate distance of 0.6a.u. between the freeform surfaces.

We want to point out that it is possible with the same algorithm to design systems with caustic surfaces between the freeform surfaces for the presented examples, which can be done by mirroring the initial map ${\mathbf{u}}^{\infty }(x,y)$ at the point of origin. This transformation preserves the property [Eq. (3)] and therefore $\nabla \mathbf{v}=0$ in Eq. (11) so that the arguments given in the previous sections hold analogously for this map. Using the mirrored intial map for the root finding and the mirrored boundary conditions, Eq. (18) leads therefore to the additional solution. We want to remark that the arguments given in the previous sections, which led to the presented design algorithm, only hold for ${\mathbf{u}}^{\infty }(x,y)$ and its mirror map, since $\nabla \mathbf{v}=0$ for both of them. This does not rule out the existence of additional solutions to the freeform design problem, such as saddle-shaped solutions with vertically or horizontally crossing rays.

The initial map ${\mathbf{u}}^{\infty }(x,y)$ can also be scaled and shifted by constants, which corresponds to a scaling of the size of the target irradiance and shifting of its position, without changing Eq. (3). These transformations are of course limited by physical boundaries, e.g., in form a total internal reflection. Additionaly, single-lens systems and two-lens systems can be calculated. Hence, the design method offers the freedom to choose between different optical configurations without a recalculation of the initial map ${\mathbf{u}}^{\infty }(x,y)$, which is one of the useful properties of the presented method.

We want to emphasize that the transformations above (scaling and shifting of the target irradiance or mirroring ${\mathbf{u}}^{\infty }(x,y)$) do not change the presented design algorithm. In our experience, the transformations do not introduce any additional numerical complications and are limited by physical constraints. Therefore, we are concentrating in the following on the case of noncrossing rays and source and target irradiances of the same size. The calculation of the initial map ${\mathbf{u}}^{\infty }(x,y)$ with our implementation of Sulman’s algorithm took 188 sec for “IAP” and 575 sec for “house” on an Intel Core i3 at $2\times 2.4$ Ghz with 16GB RAM.

After fixing the physical properties, the mapping has to be optimized according to Sections 3 and 4 to obtain the solution of the equation systems (15) and (16). To do that, one could use a starting value smaller/larger than the design goal ${\mathrm{OPL}}_{\text{red}}=-0.2$a.u. or ${\mathrm{OPL}}_{\text{red}}=0.2$a.u., respectively, to ensure the convergence to the solution of Eqs. (15) and (16). However, in practice we did not experience any significant benefit from using smaller/larger starting values for ${\mathrm{OPL}}_{\text{red}}$ than the design goal. For the tolerance of MATLAB’s fsolve(), the value of ${10}^{-6}$ was used. For both examples, the root-finding processes with 125,000 design variables took about 67 sec for “IAP” and 71 sec for “house.”

In Fig. 3 (“IAP”) and Fig. 4 (“house”), Eqs. (15) and (16) are plotted with the initial map ${\mathbf{u}}^{\infty }(x,y)$ and the optimized map $\mathbf{u}(x,y)$, respectively, to compare the quality of both mappings. Whereas the solution of the Jacobian Eq. (15) and therefore the local energy-preserving property remains nearly the same, the solution of the mapping Eq. (16) improves drastically in both cases.

 

Fig. 3. To compare the quality of the initial mapping ${\mathbf{u}}^{\infty }(x,y)$ and optimized mapping $\mathbf{u}(x,y)$, they are plugged into the Jacobian Eq. (15). The Jacobian equation for “IAP” (a) before the optimization with ${\mathbf{u}}^{\infty }(x,y)$ ($\mathrm{rms}=1.1730·{10}^{-6}$) and (b) after the optimization with $\mathbf{u}(x,y)$ ($\mathrm{rms}=1.2387·{10}^{-6}$). According to the rms values, there is a minimal decrease of quality of the local energy conservation property. The ray-mapping condition [Eq. (16)] (c) before the root finding with ${\mathbf{u}}^{\infty }(x,y)$ ($\mathrm{rms}=0.0164$) and (d) after the root finding with $\mathbf{u}(x,y)$ ($\mathrm{rms}=1.5199·{10}^{-8}$). The decrease of the rms and peak-to-valley leads to an approximate path-independent integration of the map.

Download Full Size | PDF

 

Fig. 4. To compare the quality of the initial mapping ${\mathbf{u}}^{\infty }(x,y)$ and optimized mapping $\mathbf{u}(x,y)$, they are plugged into the Jacobian Eq. (15). The Jacobian equation for “house” (a) before the optimization with ${\mathbf{u}}^{\infty }(x,y)$ ($\mathrm{rms}=5.6635·{10}^{-7}$) and (b) after the optimization with $\mathbf{u}(x,y)$ ($\mathrm{rms}=6.7198·{10}^{-7}$). There is a minimal decrease of quality of the local energy-conservation property. The ray-mapping condition [Eq. (16)] (c) before the root-finding with ${\mathbf{u}}^{\infty }(x,y)$ ($\mathrm{rms}=0.0206$) and (d) after the root finding with $\mathbf{u}(x,y)$ ($\mathrm{rms}=7.4189·{10}^{-8}$).

Download Full Size | PDF

The construction of the surfaces is done by the integration of Eqs. (9) and (13), first from $(x,y)=(-0.5,-0.5)$ along the $x$ direction and then along the $y$ direction. According to the arguments given in the previous sections, the corresponding surface will fulfill Eq. (5) and therefore be continuous. The system layouts with the freeform surfaces for both examples together with a few rays can be seen in Fig. 5.

 

Fig. 5. Layout of the lens system with freeform surfaces $z_{\mathrm{I}}(x,y)$ and $z_{\mathrm{II}}(x,y)$ mapping a collimated Gaussian beam onto a collimated beam (a) with the “IAP” irradiance distribution by a two-lens system with plane surfaces as the first and last surface and (b) with the “house” target irradiance distribution by a single optical element with two freeform surfaces.

Download Full Size | PDF

To evaluate the improvement by the optimized map, the freeform surfaces are calculated from the initial mapping ${\mathbf{u}}^{\infty }(x,y)$ and the optimized map $\mathbf{u}(x,y)$ and imported into ray-tracing software for the calculation and the comparison of the illumination patterns and the wavefronts. The results of the ray tracing with $200·{10}^6$ rays can be seen in Fig. 6 (“IAP”) and Fig. 7 (“house”). For quantitative comparisons of the results produced by ${\mathbf{u}}^{\infty }(x,y)$ and $\mathbf{u}(x,y)$, we calculated the rms, the peak-to-valley value $pv$ of the absolute difference $\mathrm{\Delta }{I}_T(x,y)≔(I_T(x,y)-I_{T,RT}(x,y))$ between predefined distribution $I_T(x,y)$ from Fig. 2 and the ray-tracing illumination patterns $I_{T,RT}(x,y)$, the correlation coefficient ${\mathrm{corr}}_{I_T}$ between $I_T(x,y)$ and $I_{T,RT}(x,y)$, and the rms values of the optical path difference. The corresponding values are listed in Tables 1 and 2. The energy efficiency for the example “IAP” was 99.9% and for “house” 99.7%.

Table 1. Comparison of $\mathsf{\Delta }{\mathsf{I}}_{\mathsf{T}}$ and OPD for Example “IAP”

View Table | View all tables in this article

Table 2. Comparison of $\mathsf{\Delta }{\mathsf{I}}_{\mathsf{T}}$ and OPD for Example “House”

View Table | View all tables in this article

 

Fig. 6. Results from the ray-tracing evaluation for the predefined irradiance distribution “IAP.” Output irradiance distribution from the ray tracing using surfaces from (a) the initial map ${\mathbf{u}}^{\infty }(x,y)$ and (b) the optimized map $\mathbf{u}(x,y)$. Absolute difference $\mathrm{\Delta }{I}_T(x,y)$ between predefined [Fig. 2(b)] and irradiance distribution from the ray tracing using surfaces from (c) the initial map ${\mathbf{u}}^{\infty }(x,y)$ and (d) the optimized map $\mathbf{u}(x,y)$. Optical path difference from the ray tracing with a reference wavelength of $\lambda =550\mathrm{nm}$ using surfaces from (e) the initial map ${\mathbf{u}}^{\infty }(x,y)$ and (f) the optimized map $\mathbf{u}(x,y)$. Since the maps are integrated first from $(-0.5,-0.5)$ along the $x$ direction and then along the $y$ direction, the deviations from the plane wavefront in (e) are consistent with Fig. 3(c).

Download Full Size | PDF

 

Fig. 7. Results from the ray-tracing evaluation for the predefined irradiance distribution “house.” Output irradiance distribution from the ray tracing using surfaces from (a) the initial map ${\mathbf{u}}^{\infty }(x,y)$ and (b) the optimized map $\mathbf{u}(x,y)$. Absolute difference $\mathrm{\Delta }{I}_T(x,y)$ between predefined [Fig. 2(c)] and irradiance distribution from the ray tracing with a reference wavelength of $\lambda =550\mathrm{nm}$ using surfaces from (c) the initial map ${\mathbf{u}}^{\infty }(x,y)$ and (d) the optimized map $\mathbf{u}(x,y)$. Optical path difference from the ray tracing using surfaces from (e) the initial map ${\mathbf{u}}^{\infty }(x,y)$ and (f) the optimized map $\mathbf{u}(x,y)$.

Download Full Size | PDF

These tables and the images in Figs. 6 and 7 demonstrate a significant improvement of quality of the illumination patterns and the wavefronts for both examples. Fractions of the rms value of $\mathrm{\Delta }{I}_T$ thereby result from ray-tracing noise, which can be estimated by averaging over several ray tracings with different seeds to lower the ${\mathrm{rms}}_{\mathrm{\Delta }{I}_T}$ values by approximately $0.1·{10}^{-6}$.

The freeform calculation with the initial maps ${\mathbf{u}}^{\infty }(x,y)$ shows strong deviations from the predefined specifications, whereas the optimized maps $\mathbf{u}(x,y)$ show high-quality illumination patterns and a wavefront uniformity far beyond the diffraction limit. Remaining deviations from the ideal wavefront (besides fundamental numerical limitations) are from Eq. (16), which is not fulfilled exactly, and therefore leads to an error accumulation along the integration path of Eqs. (9) and (13). The precision of the illumination pattern, on the other hand, is mainly limited by the precision of ${\mathbf{u}}^{\infty }(x,y)$. The main deviations from the predefined distribution $I_T(x,y)$ result from steep gradients, which can be seen especially for the example “IAP” in Fig. 6(d). This is in agreement with Fig. 3(b) and observations that were made in Ref. [17].

To demonstrate the applicability of the design algorithm to output beam sizes different from the input beam, we give the results for another example. Therefore, we choose the Gaussian input beam from the previous examples and the target irradiance distribution “house.” In contrast to the previous examples, we change the area of the target distribution to a square with a side length of 2a.u. and the reduced OPL to ${\mathrm{OPL}}_{\text{red}}=-0.8$a.u., which corresponds to $z_{\mathrm{II}}(-1,-1)-z_{\mathrm{I}}(-0.5,-0.5)=1.8598$a.u., according to Eq. (13). Since we have already calculated ${\mathbf{u}}^{\infty }(x,y)$ in the previous examples, we can simply scale the initial map by a factor of 2 and subsequently calculate $\mathbf{u}(x,y)$ from it. The calculation time of the root finding remains unchanged. The corresponding system layout can be seen in Fig. 8 and the results from the ray tracing in Fig. 9 and Table 3. Also for this example, the designed double freeform lens gives a high-quality target irradiance pattern [Fig. 9(b)] and a uniform wavefront [Fig. 9(f)], which shows the applicability of the presented algorithm to design examples with input and output beams of different size.

Table 3. Comparison of $\mathsf{\Delta }{\mathsf{I}}_{\mathsf{T}}$ and OPD for Example “House” with a Larger Output Beam

View Table | View all tables in this article

 

Fig. 8. Layout of the lens system with freeform surfaces $z_{\mathrm{I}}(x,y)$ and $z_{\mathrm{II}}(x,y)$ mapping a collimated Gaussian beam onto a collimated beam with the “house” target irradiance distribution by a two-lens system.

Download Full Size | PDF

 

Fig. 9. Results from the ray-tracing evaluation for the predefined irradiance distribution “house.” The side length of the square input beam is 1a.u. and for the square output beam 2a.u. Output irradiance distribution from the ray tracing using surfaces from (a) the initial map ${\mathbf{u}}^{\infty }(x,y)$ and (b) the optimized map $\mathbf{u}(x,y)$. Absolute difference $\mathrm{\Delta }{I}_T(x,y)$ between predefined [Fig. 2(c)] and irradiance distribution from the ray tracing with a reference wavelength of $\lambda =550\mathrm{nm}$ using surfaces from (c) the initial map ${\mathbf{u}}^{\infty }(x,y)$ and (d) the optimized map $\mathbf{u}(x,y)$. Optical path difference from the ray tracing using surfaces from (e) the initial map ${\mathbf{u}}^{\infty }(x,y)$ and (f) the optimized map $\mathbf{u}(x,y)$.

Download Full Size | PDF

6. CONCLUSION

A design method for the calculation of continuous double freeform lenses for collimated beam shaping with complex irradiance patterns beyond the paraxial approximation was presented. The method is based on the ray-mapping condition [Eq. (16)], which was derived from the law of refraction and the surface continuity condition in Section 2 and builds together with the Jacobian Eq. (15) a system of nonlinear PDEs for the unknown ray mapping $\mathbf{u}(x,y)$.

Due to the satisfaction of Eqs. (15) and (16) for infinite lens distances by the mapping from OMT with the quadratic cost function [Eq. (2)], this mapping serves as an ideal initial iterate for a root-finding approach for solving the system of Eqs. (15) and (16). As was shown by approximating Eqs. (15) and (16) by finite differences and using a standard root-finding scheme from MATLAB’s optimization toolbox, one can ensure a fast convergence to the solution of Eqs. (15) and (16). This was demonstrated by applying the presented method to several design examples with complex target distributions and validating the results by ray tracing. The double freeform surfaces showed thereby high accuracy for the irradiance patterns and the wavefront, which was assessed by the calculation of the corresponding rms values of the normalized differences.

Further improvements can be made by using OMT methods for more complex boundaries of the source and target distributions, which requires the replacement of Sulman’s method (with e.g., [29]) and the generalization of Eq. (18) to more complex boundary shapes. The scalability of the distance of the initial map ${\mathbf{u}}^{\infty }(x,y)$ by ${\mathrm{OPL}}_{\text{red}}$ to the solution of Eqs. (15) and (16) suggests the application of, e.g., the Newton algorithm for a faster root finding.

In our future research, we want to generalize the presented approach to double freeform surfaces for wavefronts different from the plane case, such as, e.g., spherical wavefronts.