1. Introduction

The problem of designing freeform optical elements generating light beams with prescribed irradiance distribution and wavefront has a wide range of applications including the design of lighting and backlight systems, laser beam processing and shaping, and optical lithography, among others. In the present work, we consider the problem of collimated beam shaping in the geometrical optics approximation. This problem consists in calculating a refracting (or reflecting) optical element, which transforms a given incident beam with a plane wavefront into an output beam with a prescribed irradiance distribution and a plane wavefront. Simultaneous control of both the irradiance distribution and the wavefront of the output beam requires the use of an optical element with two “working” optical surfaces.

In the geometrical optics approximation, the problem of collimated beam shaping can be reduced to finding a solution to a complex, fully nonlinear differential equation (NDE) of Monge–Ampère type. For the solution of this NDE, various finite-difference methods have been proposed [1–7]. The formulation of the problem of designing an optical element as an NDE assumes that the calculated surfaces of the optical element are smooth. This imposes a number of restrictions on the class of irradiance distributions that can be generated. In particular, an optical element with smooth surfaces cannot generate an irradiance distribution defined on a disconnected or a multiply connected region, or on a region with complex non-smooth boundary [8,9].

In the basic theoretical works [10–13], the problem of collimated beam shaping (PCBS) was formulated as a Monge–Kantorovich mass transportation problem (MTP) for the cases of a refractive optical element [10,12,13] and a system of two mirrors [11]. This MTP describes a problem of calculating an integrable ray mapping (i. e. a mapping relating the coordinates of the rays incident on an optical element to the coordinates of the refracted (or reflected) rays in the target plane), which provides the formation of an output beam with the required irradiance distribution. The calculation of an optical element in this MTP-based approach can be reduced to finding a solution to a linear programming problem [9] or a linear assignment problem [8]. In contrast to the NDE-based formulation, the MTP-based formulation of the PCBS enables calculating optical elements with continuous piecewise-smooth optical surfaces generating prescribed irradiance distributions defined on disconnected regions and on regions with complex non-smooth boundaries [8,9].

One of the widely used methods of designing refracting and reflecting optical surfaces is the supporting quadric method (SQM) [14–20]. It was initially proposed by V. Oliker and co-workers for solving the problem of calculating mirrors generating prescribed irradiance distributions [14–16]. The method was then generalized to the case of refracting optical surfaces [18–20]. The main advantages of the SQM are its versatility and simplicity. In the method, the prescribed irradiance distribution is approximated by a discrete distribution defined on a finite set of $N$ points. Then, the optical surface is represented as a set of $N$ segments of quadrics, each focusing the incident beam to one of these points. Depending on the problem being solved, paraboloids, ellipsoids, hyperboloids or more complex surfaces (e. g., Cartesian ovals [18]) are used as quadrics. In particular, in the problem of calculating a mirror generating a given discrete irradiance distribution in the far field, the mirror surface is represented as a set of segments of paraboloids. The calculation of the parameters of the quadrics is carried out iteratively. Note that in the above-mentioned problem of calculating a mirror, the convergence of this iterative method has been strictly proved [14].

In the present work, we consider a version of the SQM for designing a refractive optical element with two working surfaces for collimated beam shaping. The proposed method is simple to implement and is based on just a few main equations. We demonstrate that despite its simplicity, the method has firm theoretical foundations. Using the MTP-based formulation of the problem, we show that the proposed SQM is, in its essence, a gradient descent method of maximizing a concave function being a discrete analogue of the Lagrange functional in the MTP. It is important to note that any maximum (any critical point) of this function provides a solution to the PCBS. Therefore, the proposed version of the SQM does not suffer from the “false” convergence to a local maximum, which is typical for gradient methods. We illustrate the capabilities of the proposed method by designing optical elements confirming its high performance. The presented examples demonstrate the possibility of calculating optical elements with continuous piecewise-smooth surfaces generating prescribed irradiance distributions in doubly connected regions.

The paper is organized as follows. In Section 2, we formulate the problem of nonimaging optics being solved. In Section 3, we represent the surfaces of the designed element as envelopes of quadrics. Section 4 is dedicated to the description of the proposed version of the SQM. In Sections 5, we establish the connection between the SQM and the gradient method of maximizing a certain concave function associated with the MTP. In Section 6, we present two design examples illustrating the efficiency of the method. Section 7 concludes the paper.

2. Problem statement

The geometry of the considered problem is presented in Fig. 1(a). A light beam, which has a plane wavefront parallel to the plane $z=0$, impinges on an optical element with two refracting surfaces [Fig. 1(a)]. The irradiance distribution of the incident beam is defined by $I(\mathbf {x})$, $\mathbf {x} = (x_1, x_2)\in G$, where $G$ is a certain domain in the plane $z=0$. The first surface $R_1$ of the optical element is defined as $z=u_1 (\mathbf {x})$, $\mathbf {x}\in G$. It is convenient to define the second surface $R_2$ with respect to a certain “output” plane $z=f>0$, in which the irradiance distribution of the transmitted beam will be defined. Let us denote by $\mathbf {y} = (y_1, y_2)$ the Cartesian coordinates in this plane and define the second surface through the function $u_2 (\mathbf {y})$ as $z=f-u_2 (\mathbf {y})$. We assume that the refractive index of the medium between the surfaces [at $u_1 (\mathbf {x})\le z\le f-u_2 (\mathbf {x})$] equals $n_1$, and the refractive index of the surrounding medium [at $z < u_1 (\mathbf {x})$ and at $f-u_2 (\mathbf {x}) < z$] equals $n_0$.

 

Fig. 1. (a) Geometry of the problem of collimated beam shaping. (b) Representation of the first surface as the envelope of a family of surfaces focusing light to the points of the second surface (two-dimensional case). The envelope is shown with a dashed black curve. Solid colored curves show the cross-sections of different surfaces of the family.

Download Full Size | PDF

An incident ray from a point $\mathbf {x}\in G$ of the plane $z=0$ is successively refracted by the surfaces $R_1$ and $R_2$ and arrives at some point $\mathbf {y}$ of the “target” domain $D$ in the output plane $z=f$ [Fig. 1(a)]. In this sense, the surfaces $R_1$ and $R_2$ define a ray mapping $\gamma : G\to D$, which relates the coordinates of the refracted rays $\mathbf {y}=\gamma (\mathbf {x})$ to the coordinates of the incident rays. This mapping determines the irradiance distribution $L(\mathbf {y})$, $\mathbf {y}\in D$ generated in the plane $z=f$ according to the light flux conservation law:

The problem of collimated beam shaping, which we consider in the present work, is formulated in the following way. Assuming that the irradiance distribution of the incident beam $I(\mathbf {x})$, $\mathbf {x}\in G$ and the required irradiance distribution $L(\mathbf {y})$, $\mathbf {y}\in D$ in the output plane $z=f$ [Fig. 1(a)] are given, it is necessary to find such functions $u_1(\mathbf {x})$ and $u_2 (\mathbf {y})$ so that the transmitted light beam would have the given irradiance distribution $L(\mathbf {y})$ and a plane wavefront parallel to the plane $z=f$.

3. Representation of the optical surfaces

Let us consider a ray from a point $\mathbf {x}$ of the plane $z=0$, which, after being refracted by the surfaces $R_1$ and $R_2$, arrives at a point $\mathbf {y}=\gamma (\mathbf {x})$ of the plane $z=f$. Since the wavefronts of both the incident and the output beams are plane, the optical path length $F$ between these points is constant and does not depend on $\mathbf {x}$. This condition can be written as

Let us show that the surfaces $R_1$ and $R_2$ of the optical element can be represented as envelope surfaces of a special kind. Indeed, for differentiable surfaces $R_1$ and $R_2$, the ray mapping $\mathbf {y}=\gamma (\mathbf {x})$ can be described in the following way. The ray coming from a point $\mathbf {x}\in G$, after being refracted at the point $\left (\mathbf {x}, u_1(\mathbf {x})\right )$ of the surface $R_1$, goes to the point $\left (\mathbf {y}, f-u_2 (\mathbf {y})\right )$ of the surface $R_2$. Therefore, the surface $R_1$ can be represented as the envelope of a family of surfaces, each focusing the incident plane beam to one of the points $\left (\mathbf {y}, f-u_2 (\mathbf {y})\right )$ of the surface $R_2$ [Fig. 1(b)] [9].

The equation of the refracting surface $z=u_1^{\mathbf {y}} (\mathbf {x})$ focusing the incident beam to a certain point $\left (\mathbf {y}, f-u_2 (\mathbf {y})\right )$ can be obtained by requiring the optical path length from the points $\mathbf {x}$ of the plane $z=0$ to the point $\left (\mathbf {y}, f-u_2 (\mathbf {y})\right )$ to be constant and equal to $F-n_0 u_2 (\mathbf {y})$. Note that Eqs. (4) and (5) describe exactly this condition, and therefore the function $u_1^{\mathbf {y}} (\mathbf {x})$ has the form

Considering $\mathbf {y} = (y_1, y_2)$ as parameters, we can treat the equation $z=u_1^{\mathbf {y}} (\mathbf {x})$ as a two-parameter family of surfaces. The first surface of the designed optical element, which is described by the function $u_1 (\mathbf {x})$, is the envelope of this family [9] and is defined by the following equations [22]:

Similarly, the second surface $R_2$ can be considered as the envelope of a family of surfaces

In what follows, we will consider a particular case of the envelope of Eqs. (7a) and (7b), when the condition of Eq. (7b) corresponds not just to some critical point of the function ${\mathcal {K}}(\mathbf {x}, \mathbf {y})-u_2 (\mathbf {y})$ with respect to $\mathbf {y}$, but to its minimum. Such a choice provides the so-called Galilean configuration of rays [10]. In this case, the envelope surface defined by Eqs. (7a) and (7b) takes the form

Equations (10) and (11) make it possible to formulate the problem of collimated beam shaping as a problem of calculating such a function $u_2 (\mathbf {y})$, so that the corresponding ray mapping (11) satisfies the light flux conservation law of Eq. (1) or Eq. (2). Let us note that not every mapping satisfying the light flux conservation law is integrable, i. e. can be implemented by a pair of optical surfaces. In order for the mapping $\gamma$ to be integrable, it is necessary for the function $u_2 (\mathbf {y})$ defining it to exist and to satisfy Eq. (7b). Therefore, below we will refer to Eq. (7b) as the integrability condition.

4. Supporting quadric method

The supporting quadric method (SQM) [14–20] is widely used for designing mirrors and refracting surfaces generating prescribed irradiance distributions. In the method, the required continuous irradiance distribution is approximated by a discrete distribution consisting of a finite set of points. The designed optical surface is assumed to be piecewise-smooth and consisting of segments of quadrics. The parameters of the quadrics are found iteratively from the condition of obtaining the required light flux values at the points of the discrete distribution.

In this section, we consider a version of the SQM for solving the problem of collimated beam shaping [Fig. 1(a)]. Strictly speaking, a rigorous discussion of the SQM requires the introduction of the so-called “weak solution” concept, which makes it possible to correctly define the solution in the cases, in which the ray mapping is discontinuous on a set of measure zero [9]. However, for the sake of simplicity, we will concentrate on the practical side of the method. Let the prescribed continuous irradiance distribution $L(\mathbf {y})\in D$ be approximated by a discrete distribution $L_{\textrm {d}} (\mathbf {y})=\sum _{i=1}^{N} L_i \delta (\mathbf {y}-\mathbf {y}_i)$ (a sum of Dirac delta functions) concentrated at $N$ points $\mathbf {y}_i$, $i=1, \ldots , N$. In this case, the function $u_2 (\mathbf {y})$ in Eqs. (10) and (11) is replaced by a set of $N$ numbers $u_{2, i}$, $i=1, \ldots , N$. It is these values that are varied in the method. In this case, refractive surface $R_1$ is a piecewise-smooth surface consisting of segments of refractive lenses of Eq. (6) (ellipsoids or hyperboloids) focusing the incident plane beam to the points $(\mathbf {y}_i, f-u_{2, i})$, $i=1, \ldots , N$. This surface is schematically shown in Fig. 1(b) with bold segments of the colored curves showing the ellipsoids focusing the incident light to the individual points of the second optical surface. According to Eq. (10), the function $u_1 (\mathbf {x})$ defining the surface $R_1$ has the form

In practical calculations, the integral in the right-hand side of Eq. (15) can be calculated numerically using the ray tracing method, which is schematically shown in Fig. 1(b). In this approach, the incident beam is approximated by a set of $N_\textrm {tr}$ rays going from the nodes of a certain grid $\mathbf {x}_j$, $j=1, \ldots , N_\textrm {tr}$ in the plane $z=0$. These rays “carry” the light flux values $I_j$, $j=1, \ldots , N_\textrm {tr}$, which are proportional to the values of the irradiance of the incident beam at the grid nodes. According to Eq. (13), for each ray coming from a certain point $\mathbf {x}_j$, we find an index $i(j)=\mathop {\arg \min }\nolimits _i \left [{\mathcal {K}}(\mathbf {x}_j, \mathbf {y}_i)-u_{2,i} \right ]$ corresponding to the lens, which the ray intersects first [Fig. 1(b)]. This ray is directed to the corresponding point $\left (\mathbf {y}_{i(j)}, f-u_{2,i(j)} \right )$. Then, at the points of the second optical surface, the light flux values are calculated, which are equal to the sums of the light flux values of the rays focused to these points. In Fig. 1(b), these values are shown by numbers near each of the points of the “discretized” second optical surface. It is easy to see that with an increase in the $u_{2, i}$ value (in the case of fixed $u_{2, j}, j\ne i$), the domain $C(\mathbf {y}_i)$ increases, so that the light flux $L_{\textrm {est}, i}$ increases as well. Thus, if at some point $u_{2, i}$, the required light flux is greater than the obtained one [$L_{\textrm {d}} (\mathbf {y}_i)=L_i >L_{\textrm {est}, i}$], it is necessary to increase the $u_{2i}$ value (and to decrease it in the opposite case). It is natural to choose the magnitude of this increase or decrease proportional to the difference between the required and calculated light flux values $\left (L_i -L_{\textrm {est}, i} \right )$ [20]. In this way, the values $u_{2, i}$ can be found iteratively. At each iteration, the generated distribution $L_{\textrm {est}, i}$, $i=1, \ldots , N$ is calculated, and the $u_{2, i}$ values are updated according to the rule

It is important to note that according to the light beam reversibility principle, the considered PCBS can be formulated in a “reverse” form. In this case, the propagation direction is reversed and the input and required irradiance distributions are interchanged. In such a reverse formulation, one can discretize the irradiance distribution $I(\mathbf {x})$ [instead of $L(\mathbf {y})$] and then use an equation similar to Eq. (16) for the calculation of the values $u_{1,i}$ representing the first surface.

5. Connection of the SQM with a gradient method

5.1 Mass transportation problem

In Ref. [10], it was shown that in the considered collimated beam shaping problem, the calculation of a ray mapping satisfying the light flux conservation law (1) can be reduced to finding a solution to a Monge–Kantorovich mass transportation problem (MTP) with the cost function ${\mathcal {K}}(\mathbf {x}, \mathbf {y})$ of Eq. (5). The strict formulation of this fact is given by the following theorem:

Theorem 1. Let the mapping $\gamma : G\to D$ provide a minimum to the functional

Proof. Let us consider the Lagrange functional for the considered variational problem [25]:

The variations of the functional ${\mathcal {H}}(T, \lambda )$ can be easily calculated. From Eq. (20), one can easily show that the condition $\delta _{\lambda } {\mathcal {H}}=0$ of zero variation with respect to $\lambda$ is equivalent to the light flux conservation law [26]. To calculate the variation $\delta _{T} {\mathcal {H}}$ with respect to $T$, let us split the integral in the right-hand side of Eq. (20) into two integrals and make the change of variables $\mathbf {y}=T(\mathbf {x})$ in the second integral:

The obtained equations are equivalent to the integrability condition (7b). According to the conditional extremum theorem [25], if the mapping $\gamma$ provides a conditional extremum to the problem of Eqs. (17) and (18), then there exists a function $\lambda (\mathbf {y})$, which satisfies Eqs. (22). This means that the mapping $\gamma$ is integrable, i. e. can be implemented by a system of two refracting surfaces described by the functions $u_1 (\mathbf {x})$ and $u_2 (\mathbf {y})$, where $u_2 (\mathbf {y})\equiv \lambda (\mathbf {y})$, and the function $u_1 (\mathbf {x})$ is defined through $u_2 (\mathbf {y})$ using Eq. (10). This concludes the proof of the theorem.

From the presented proof, it is evident that the solution of the collimated beam shaping problem can be obtained from any critical point of the Lagrange functional (20), and not only from the point corresponding to the minimum of the functional (17). In particular, the mapping $\gamma$, which provides a maximum to the functional (17), also gives a solution of the problem. From the form of the function ${\mathcal {K}}(\mathbf {x}, \mathbf {y})$ [see Eq. (5)], it is clear that the conditions of minimization and maximization of the functional (17) correspond to the Galilean and Keplerian configurations of rays, respectively [8,9].

Let us note that the ray mapping appears under the integral in the Lagrange functional (20). This makes it possible to use the MTP approach for calculating refracting surfaces implementing not only smooth mappings $\gamma$, but also mappings, which are discontinuous on a set of measure zero. In the case of a discontinuous mapping $\gamma$, the optical surface $u_1 (\mathbf {x})$ calculated using Eq. (10) will be piecewise-smooth.

5.2 Modification of the Lagrange functional

From the proof of Theorem 1, it follows that if a pair $(T, \lambda )$ is a critical point of the Lagrange functional ${\mathcal {H}} (T, \lambda )$, then the mapping $\mathbf {y}=T(\mathbf {x})$ satisfies the light flux conservation law and is defined by a system of refracting surfaces $u_1 (\mathbf {x})$ and $u_2 (\mathbf {y})$, where $u_2 (\mathbf {y})=\lambda (\mathbf {y})$, and $u_1 (\mathbf {x})$ is found from Eq. (10). At a critical point, the mapping $\mathbf {y}=T(\mathbf {x})$ and the function $\lambda (\mathbf {y})$ are related by the integrability condition (22). Therefore, the mapping can be expressed through the function $\lambda (\mathbf {y})$ using Eq. (22). Similarly to the discussion above, where we restricted our consideration to the envelopes defined by Eq. (10), in what follows, we consider the critical points of a special kind. Namely, we assume that the condition of the equality of the derivatives to zero in Eq. (22) corresponds not just to some critical point of the function ${\mathcal {K}}(\mathbf {x}, \mathbf {y})-\lambda (\mathbf {y})$, but to a minimum point. In this case, the mapping $T$ can be written in a form similar to the expression for the ray mapping (11):

5.3 Generalized Voronoi cells and gradient descent method

In this section, we establish the connection between the functional (24) and the supporting quadric method (16). For this, let us consider the form of the functional (24) when the continuous irradiance distribution $L(\mathbf {y})$ is approximated by a discrete distribution $L_{\textrm {d}} (\mathbf {y})=\sum _{i=1}^{N} L_i \delta \left (\mathbf {y}-\mathbf {y}_i \right )$ defined on $N$ points $\mathbf {y}_i$, $i=1, \ldots , N$. In this case, the function $\lambda (\mathbf {y})$ is replaced by a vector $\boldsymbol \lambda = (\lambda _1, \ldots , \lambda _N)$, where $\lambda _i =\lambda (\mathbf {y}_i)$, and the mapping $T^{\lambda }$ becomes a mapping of the domain $G$ to a discrete set of points $\mathbf {y}_i$, $i=1, \ldots , N$.

Definition. Let us define the generalized weighted Voronoi cells as the following domains $C^{\boldsymbol \lambda } (\mathbf {y}_i)$: a point $\mathbf {x}\in G$ belongs to the domain $C^{\boldsymbol \lambda } (\mathbf {y}_i)$ if and only if

 

Fig. 2. Partition (a) of the domain $G$ into equal-area generalized Voronoi cells $C^{\boldsymbol \lambda } (\mathbf {y}_i)$ associated with the points $\mathbf {y}_i$ (b).

Download Full Size | PDF

In the case of a discrete distribution $L_{\textrm {d}} (\mathbf {y})$, the functional ${\mathcal {H}}(\lambda )$ becomes a function defined on ${\mathbb {R}}^{N}$:

In the Appendix, we consider the calculation of the partial derivatives of the function (26) for an arbitrary function ${\mathcal {K}}(\mathbf {x}, \mathbf {y})$. According to the results presented in the Appendix, these derivatives have the form

Note that since the function $\lambda (\mathbf {y})$ coincides with the function $u_2 (\mathbf {y})$ defining the surface $R_2$, then $u_{2,i}=\lambda _i$ and the generalized Voronoi cells $C^{\boldsymbol \lambda } (\mathbf {y}_i)$ coincide with the domains $C (\mathbf {y}_i)$ [Eq. (14)], on which the SQM quadrics are defined, which focus the incident beam to the points $\mathbf {y}_i$, $i=1, \ldots , N$. The values $L_{\textrm {est}, i} = \int _{C^{\boldsymbol \lambda } (\mathbf {y}_i)} I(\mathbf {x})\mathrm {d}{\mathbf {x}}$ are equal to the light flux focused to the points $\mathbf {y}_i$, $i=1, \ldots , N$. According to Eq. (27), at a critical point of the function $h(\boldsymbol \lambda )$, the following conditions hold:

It is important to note that the above theoretical analysis is different from the theoretical results of papers [10,12], which, in the discrete case, allow one to reduce the PCBS to a linear programming problem (LPP). In contrast, in the present work we reduce the PCBS to an unconstrained optimization problem, which consists in maximizing the concave nonlinear function $h(\boldsymbol \lambda )$ being a discrete analogue of the modified Lagrange functional (24). This approach, in comparison with the LPP-based formulations in Refs. [10,12], turns out to be much simpler in terms of the computational complexity. Indeed, the LPP representing the PCBS has large dimensions and, in most cases, cannot be solved using the commonly available LPP solvers (e.g. SOPLEX, MOSEK, or linprog MATLAB routine). For example, if both surfaces of the optical element are defined on grids with the dimensions of $100\times 100$ containing $N=10000$ points each (in this case, the incident beam and the generated collimated beam are also defined on grids with the same dimensions), the number of constraints in the mentioned LPP will be equal to $N^{2}=10000^{2}=10^{8}$ [10]. The LPP with such a large number of constrains is extremely complex for the existing LPP solvers. Thus, the “direct” application of the LPP-based technique of Refs. [10,12] is possible only for rather coarse meshes (with the sizes of about $50\times 50$), which is insufficient for many practical problems. The design examples presented in the next section demonstrate that the proposed SQM version of Eq. (16) [the gradient descent method of maximizing the function (26)] enables designing high-performance optical elements generating collimated beams with required irradiance distributions defined on the grids with sizes $100\times 100$ and $150\times 150$.

6. Design examples

6.1 Generation of a “square-frame” beam with uniform irradiance

As a first example, let us consider the calculation of an optical element transforming a beam with a circular cross-section and uniform irradiance

For the calculation of the optical surfaces, the required irradiance distribution of the output beam (30) was approximated by a discrete distribution defined on a square grid with spacing $\Delta _y = 0.0173\textrm {~mm}$. In this case, the discrete distribution contains $N=100^{2}$ points. For the calculation of the second surface (of the values $u_{2,i}$, $i=1,\ldots ,N$), we used the SQM version described in Section 4. Under this approach, at each iteration the light flux at the points $\left (\mathbf {y}_i, f-u_{2,i} \right )$, $i=1,\ldots ,N$ is calculated using the ray tracing technique, and then the values $u_{2,i}$ are updated using Eq. (16). In the ray tracing method, the incident beam was approximated by a set of $N_\textrm {tr} =2N=20000$ rays defined on a square grid in the plane $z=0$. In the considered example, 25000 iterations were performed. The total calculation time on a desktop computer (Intel Core i9-7940X CPU) was of about 6 minutes.

After the calculations, the second surface defined by a discrete set of values $u_{2,i}$, $i=1,\ldots ,N$ was approximated using the least-squares approach by a system of two-dimensional B-splines [30]. Then, from the calculated second surface, the first surface was reconstructed using Eq. (10). The obtained optical surfaces are shown in Fig. 3(a). It is interesting to note that in the vicinity of the origin the first surface has sharp bends (derivative discontinuities) and resembles a tetrahedral pyramid. This is due to the fact that the central part of the first surface maps the central part of the incident beam (the neighborhood of the point $\mathbf {x}=0$) to the inner contour of the square-frame. This mapping is discontinuous. The possibility to design piecewise-smooth surfaces implementing discontinuous ray mappings is an important advantage of the SQM. In this case, the existing methods based on finding a numerical solution to a partial differential equation of elliptic type [1–7] are either inapplicable or numerically unstable.

 

Fig. 3. Optical surfaces of the designed element (a) and generated normalized irradiance distributions in the planes $z=10\textrm {~mm}$ (b) and $z=30\textrm {~mm}$ (c). Red and blue lines show the central irradiance cross-sections.

Download Full Size | PDF

Figures 3(b) and 3(c) show the calculated normalized irradiance distributions generated by the designed optical element in the output plane $z=f=10\textrm {~mm}$ and in the plane $z=30\textrm {~mm}$. Note that these distributions were calculated on a denser grid with the spacing $\Delta _{y} /4=0.0043\textrm {~mm}$ by tracing $10^{7}$ rays. The presented simulation results demonstrate the formation of a square-frame beam with the required dimensions and almost uniform irradiance. The normalized root-mean-square deviations of the obtained irradiance distributions from the uniform one amount to 3.3% and 3.8% at $z=10\textrm {~mm}$ and $z=30\textrm {~mm}$, respectively. Let us note that the dimensions of the square-frame beam at different distances from the optical element remain the same, therefore, the wavefront of the generated output beam is plane. Indeed, the calculated root-mean-square (RMS) deviation of the optical path length from the constant design value $F=11\textrm {~mm}$ amounts to $9\textrm {~nm}$ at $z=10\textrm {~mm}$.

6.2 Generation of an annular beam with non-uniform irradiance

As the second example, we designed an optical element transforming the incident beam (29) with the radius $R=1\textrm {~mm}$ into an annular beam with a non-uniform irradiance distribution

The calculated surfaces of the optical element are shown in Fig. 4(a). As in the first example, the first surface has sharp bends near the origin of coordinates, which is due to the discontinuity of the ray mapping, which maps the neighborhood of the point $\mathbf {x}=0$ to the inner boundary of the ring. Figures 4(b) and 4(c) show the simulated normalized irradiance distributions generated by the designed element in the output plane $z=f=10\textrm {~mm}$ and in the plane $z=30\textrm {~mm}$. As in the previous example, the distributions were calculated on a grid with $0.0043\textrm {~mm}$ spacing by tracing $10^{7}$ rays. The presented results demonstrate the formation of an annular beam with the required dimensions. The normalized root-mean-square deviations of the calculated irradiance distributions from the required distribution defined by Eq. (31) amount to 6.2% and 6.7% at $z=10\textrm {~mm}$ and $z=30\textrm {~mm}$, respectively. The wavefront of the generated beam is almost plane and the calculated root-mean-square deviation of the optical path length from the constant design value $F=11\textrm {~mm}$ amounts to $13\textrm {~nm}$ at $z=10\textrm {~mm}$.

 

Fig. 4. Optical surfaces of the designed element (a) and generated normalized irradiance distributions in the planes $z=10\textrm {~mm}$ (b) and $z=30\textrm {~mm}$ (c). Red and blue lines show the central irradiance cross-sections.

Download Full Size | PDF

6.3 Discussion of the convergence criteria

The SQM convergence is usually controlled by using criteria that represent the difference between the calculated and prescribed discrete distributions (e. g., the RMS deviation) [15,20]. However, such a choice is not optimal in the problem under consideration. Indeed, according to the results presented in Section 5.3, the SQM of Eq. (16) is a gradient method of maximizing a concave function $h(\boldsymbol \lambda )$ [Eq. (26)], which is a discrete counterpart of the Lagrange functional. To illustrate the concavity of this function, we chose two random vectors, $\boldsymbol \lambda _1$ and $\boldsymbol \lambda _2$, and plotted the function along a line segment connecting them: $h(t \boldsymbol \lambda _1 + (1-t)\boldsymbol \lambda _2)$, $t\in [0, 1]$ [see the upper panel in Fig. 5(a)]. For comparison, the lower panel of Fig. 5(a) shows the RMS deviation between the vectors consisting of the prescribed, $\mathbf {L}$, and the calculated, $\mathbf {L}_\textrm {est}(t \boldsymbol \lambda _1 + (1-t)\boldsymbol \lambda _2)$, light flux values. The presence of several local minima and noisy behavior of the RMS deviation make it a poorer criterion for controlling the convergence of the iterative process than $h(\boldsymbol \lambda )$ or a certain monotonic function of it.

 

Fig. 5. (a) Convergence criteria calculated at $\boldsymbol \lambda = t \boldsymbol \lambda _1 + (1-t)\boldsymbol \lambda _2$, $t\in [0, 1]$: $h(\boldsymbol \lambda )$ (upper panel) and RMS of the difference between the prescribed and calculated distributions (lower panel). (b) Convergence rate of the method: the difference $h_\textrm {opt} - h(\boldsymbol \lambda )$ (upper panel) and the RMS of the difference between the prescribed and calculated distributions (lower panel) vs. the iteration number.

Download Full Size | PDF

Let us illustrate this fact using as an example the optical element generating a square-frame beam (Fig. 3). The upper panel of Fig. 5(b) shows the plot of the quantity $h_\textrm {opt} - h(\boldsymbol \lambda )$, where $h(\boldsymbol \lambda )$ is calculated using Eq. (26) at each SQM iteration at $\boldsymbol \lambda =(u_{2,1}, \ldots ,u_{2,N})$, and $h_\textrm {opt}$ is the final value of this function obtained at the last iteration. From Fig. 5(b), it is evident that the difference $h_\textrm {opt} - h(\boldsymbol \lambda )$ monotonically decreases. An abrupt decrease in this value at iteration 5001 is caused by a decrease in the step $\alpha$ in Eq. (16): the first 5000 iterations were performed with the step $\alpha =1$, and then the step was decreased to $\alpha = 0.1$. A similar convergence plot for the RMS criterion is shown in the lower panel of Fig. 5(b). In contrast to $h(\boldsymbol \lambda )$, the RMS deviation graph exhibits a counterintuitive increase starting from the 5001st iteration, and is much more noisy. Therefore, the RMS criterion is not suitable for controlling the convergence, and it is the $h(\boldsymbol \lambda )$ function that should be used for this purpose.

7. Conclusion

We proposed a version of the supporting quadric method for designing refractive optical elements for collimated beam shaping. We demonstrated that this method is, in its essence, the gradient method of maximizing a concave function, which is a discrete analogue of the Lagrange functional in a corresponding mass transportation problem. We established that any maximum of this function provides a solution to the problem of collimated beam shaping. Therefore, the proposed method does not suffer from the convergence to a local extremum, which is typical for the “conventional” gradient methods. Another important advantage of the method is simplicity, since its implementation is based on only two simple equations: Eqs. (13) and (16). It is important to note that instead of the simple gradient descent method of Eq. (16), one can use more sophisticated convex optimization algorithms (e. g. Nesterov’s accelerated gradient method) to obtain a higher convergence rate.

As examples, we designed refractive optical elements, which transform a circular incident beam into a square-frame beam with uniform illuminance and to an annular beam with non-uniform illuminance. The presented examples demonstrate the possibility of designing optical elements with continuous piecewise-smooth surfaces and show high performance of the proposed method.

The authors believe that the connection between the SQM and the gradient method of maximizing a function related to a certain MTP can also be obtained for other problems of nonimaging optics, which can be formulated as an MTP. In particular, similar versions of the SQM can be formulated for the problems of calculating mirrors and refracting surfaces generating prescribed irradiance distributions in the far field. The differences will only appear when defining the generalized weighted Voronoi cells, the form of which is determined by the cost function, depending on the problem of nonimaging optics being solved.