1. Introduction

The problem of calculating optical elements with freeform surfaces that generate light beams with a prescribed irradiance distribution and a prescribed wavefront has many practical applications, including the design of lighting and backlight systems, laser beam shaping, in particular, in lithography systems, laser printing, optical memory, etc. [1,2]. Simultaneous control of both the irradiance distribution and the wavefront of the output beam requires the use of an optical element with at least two “working” surfaces (refracting or reflecting).

In the geometrical optics approximation, the problem of calculating an optical element that generates prescribed irradiance distribution and wavefront can be reduced to finding a solution to a nonlinear differential equation (NDE) of elliptic type [3–10]. Despite the fact that various finite-difference methods have been proposed to solve such NDEs [3–9], the calculation of optical elements within this approach is complex and has limitations. In particular, the formulation of the problem of calculating an optical element in the form of an NDE assumes that the calculated surfaces of the optical element are smooth. The requirement for smoothness significantly limits the class of irradiance distributions that can be generated by the optical element. For example, an optical element with smooth surfaces does not enable the formation of an irradiance distribution defined on a disconnected region, or on a region with non-smooth boundaries [11,12].

In the particular case of the so-called collimated beam shaping problem, when the wavefronts of the incident and output beams are plane, the problem of calculating the optical element can be formulated as a variational problem or, in particular, as a Monge–Kantorovich mass transportation problem (MTP) [11–14]. Such an MTP describes the calculation of an integrable ray mapping (i. e., a mapping that connects the coordinates of the rays incident on the optical element and the coordinates of the refracted (reflected) rays in the output plane), which ensures the formation of an output beam with a prescribed irradiance distribution. Calculation of an optical element within the framework of this “MTP-formulation” can be reduced to solving a linear programming problem [12,14] or a linear assignment problem [11]. The MTP-formulation of the problem, in contrast to the problem of solving an NDE, makes it possible to calculate optical elements with continuous piecewise-smooth surfaces. This enables generating beams with irradiance distributions specified in disconnected regions and in regions with complex and non-smooth boundaries [11,12]. Nevertheless, when the wavefront of the input or output beam is not plane, the problem of calculating an optical element can no longer be reduced to an MTP, so that the methods of [11,13] become inapplicable. The variational methods developed in [12,14] also assume that the incident and output beams have plane wavefronts and also become inapplicable when the output beam has a non-planar wavefront.

One of the methods widely used for calculating reflecting and refracting surfaces that generate required irradiance distributions is the supporting quadric method (SQM) [15–21]. The main advantages of the SQM are its versatility and simplicity. In this method, the required irradiance distribution is approximated by a discrete distribution defined on a finite number of $N$ points. Then, the optical surface is represented as a set of $N$ fragments of quadrics that focus the incident beam to the required points. Depending on the problem being solved, paraboloids, ellipsoids, hyperboloids, or more complex surfaces, for example, Cartesian ovals [19], are used as quadrics. In particular, in the problem of calculating a mirror generating a prescribed discrete irradiance distribution in the near field, the surface of the mirror is represented as a set of segments of ellipsoids, one focus of which coincides with the point light source, and the other with one of the points of the generated discrete distribution. The calculation of the parameters of the quadrics is carried out by an iterative method. In the above-mentioned problem of calculating the mirror, the convergence of this method has been rigorously proved [16].

In a recent work of the present authors, a version of the SQM was considered for solving the problem of collimated beam shaping, i. e. for the case of plane wavefronts of the incident and output beams [22]. It was rigorously shown that the version of the SQM proposed in [22] is, in fact, a gradient descent method for maximizing a concave function, which is a discrete analogue of the Lagrange functional in the mass transportation problem. In the present work, we consider the SQM for a more general case, when the wavefront of the output beam is not plane and is specified through the eikonal function in the output plane. Similarly to [22], the proposed method is simple to implement and is based on only a few main equations. The proposed method is illustrated by examples of the design of optical elements demonstrating its high performance. The presented examples show the possibility of calculating optical elements that generate output beams with spherical wavefronts (converging and diverging) and with prescribed irradiance distributions, including an irradiance distribution specified in a doubly-connected region in the form of a square frame. To illustrate the broad capabilities of the method, we also consider the design of an optical element that forms a light beam with a uniform irradiance distribution in the square-frame region and with a complex eikonal distribution, which provides the subsequent focusing of the output beam into a circle.

2. Problem statement

Let us consider a three-dimensional space ${\mathbb {E}}^{3}$ with the coordinates $(x_1 ,x_2 ,z)$. Let a light beam with irradiance distribution $I(\mathbf {x}),\, \, \mathbf {x}=(x_1 ,x_2)\in G$ and a plane wavefront parallel to the plane $z=0$ impinge on an optical element with two refracting surfaces [Fig. 1(a)]. The first surface $R_1$ is defined by the function $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$ located after the element. In this plane, the required irradiance distribution and the eikonal function of the transmitted optical beam will be defined [Fig. 1(a)]. Let us denote by $\mathbf {y}=(y_1 ,y_2)$ the Cartesian coordinates in the plane $z=f$. The eikonal function $\Psi (\mathbf {y}),\, \mathbf {y}\in D$ describes the required wavefront of the output beam and defines the ray directions in the plane $z=f$ as

 

Fig. 1. (a) Geometry of the problem of generating a beam with prescribed irradiance distribution and wavefront (eikonal function). For example, the case of generating a beam with a converging spherical wavefront is shown. (b) Representation of the first surface (two-dimensional cross-section shown) as an envelope of a family of surfaces (ellipsoids) focusing the incident beam to the points of the second surface. Colored lines show the two-dimensional cross-sections of several surfaces of the family.

Download Full Size | PDF

Since the ray propagation directions in the plane $z=f$ have to coincide with the directions of the rays refracted by the second surface $R_2$, let us define the second surface through the function $l(\mathbf {y})$ equal to the distance along the ray direction from the second surface to the plane $z=f$ [Fig. 1(a)]. In this case, the second surface can be represented in the following parametric form:

A ray coming out of a point $\mathbf {x}\in G$ of the plane $z=0$ is sequentially refracted by the surfaces $R_1$ and $R_2$, and arrives to a certain point $\mathbf {y}$ of the output plane $z=f$ [Fig. 1(a)]. Thus, the surfaces $R_1$ and $R_2$ define a ray mapping $\mathbf {y}=\gamma (\mathbf {x})$, which relates the coordinates of the refracted rays in the output plane to the coordinates of the incident rays. This mapping determines the irradiance distribution $L(\mathbf {y})$ generated in the plane $z=f$. This distribution can be found using the light flux conservation law:

The problem of generating prescribed irradiance distribution and wavefront considered in the present work is formulated in the following way. Given an irradiance distribution $I(\mathbf {x})$, $\mathbf {x}\in G$ of the incident beam with a plane wavefront, it is necessary to find such functions $u_1 (\mathbf {x})$ and $l(\mathbf {y})$ defining the surfaces of the optical element, so that the transmitted light beam would have certain required irradiance distribution $L(\mathbf {y}),\, \, \mathbf {y}\in D$ and eikonal $\Psi (\mathbf {y}),\, \, \mathbf {y}\in D$ in the output plane $z=f$ [Fig. 1(a)].

3. Representation of the optical surfaces

The second surface $R_2$ of the optical element is defined by Eq. (2). This equation is written taking into account the required eikonal function $\Psi (\mathbf {y})$ in the output plane $z=f$. The first surface of the element $R_1$ can be represented as an envelope of a family of surfaces of a special type. Indeed, for differentiable surfaces $R_1$ and $R_2$, the ray mapping $\mathbf {y}=\gamma (\mathbf {x})$ can be described in the following way. A ray coming out of a certain point $\mathbf {x}\in G$, after being refracted at the point $\left (\mathbf {x},u_1 (\mathbf {x})\right )$ of the surface $R_1$, arrives to the point $\vec {u}_2 (\mathbf {y})=\left (\mathbf {y},f\right )-\vec {p}(\mathbf {y})l(\mathbf {y})$ of the surface $R_2$. Therefore, the surface $R_1$ can be represented as an envelope of a family of surfaces focusing an incident plane beam to the points of the surface $R_2$ [Fig. 1(b)] [8].

The equation of a refractive surface (lens) $z=u_\textrm {lens} (\mathbf {x})$ focusing the incident beam to the point $\vec {u}_2 (\mathbf {y})=\left (u_{2,x_1 } (\mathbf {y}),\, u_{2,x_2 } (\mathbf {y}),u_{2,z} (\mathbf {y})\right )$ can be easily obtained from the condition of equality of the optical path length of the rays from the points of the plane $z=0$ to the point $\vec {u}_2 (\mathbf {y})$ to the value of the eikonal function at this point $\Psi _2 (\mathbf {y})=\Psi (\mathbf {y})-n_0 l(\mathbf {y})$. This condition has the form

By definition, the envelope surface is tangent to each of the quadrics of the family defined by Eq. (6) at a certain point. According to Eq. (7), at the point of tangency, the partial derivatives of the envelope surface are equal to the partial derivatives of the quadric. Since the directions of the refracted rays are determined by these partial derivatives, the rays refracted by the envelope surface will be directed to the points of the second surface.

Note that at a fixed value of $\mathbf {x}=\left (x_1 ,x_2 \right )$, Eq. (7) describes a critical point of the function $u_\textrm {lens} (\mathbf {x};\mathbf {y},l(\mathbf {y}))$ with respect to the variables $y_1$ and $y_2$. In what follows, we will consider a particular case of the envelope (6), (7), when the condition (7) corresponds not just to some critical point of the function $u_\textrm {lens} (\mathbf {x};\mathbf {y},l(\mathbf {y}))$, but to a minimum point (the case of a maximum point can be considered similarly). In this case, the envelope surface takes the form

Equation (8) enables expressing the ray mapping as

Since the function $u_\textrm {lens} (\mathbf {x};\mathbf {y},l(\mathbf {y}))$ depends on the function $l(\mathbf {y})$ defining the second surface, Eqs. (8) and (9) enable formulating the problem of generating prescribed irradiance distribution and wavefront as a problem of calculating such a function $l(\mathbf {y})$, so that the corresponding ray mapping (9) satisfies the light flux conservation law defined by Eq. (3) or (4).

Note that not all eikonal and irradiance distributions $\Psi (\mathbf {y})$ and $L(\mathbf {y})$ can be implemented by an optical element. This is due to the limited deviation angles of the rays (i. e. the angles between the rays incident on the surface and the corresponding refracted rays) at the surfaces of the optical element. In particular, the expression under the square root in Eq. (6) specifying the lens could become negative at large offset values $\left \| \mathbf {u}_2 (\mathbf {y})-\mathbf {x}\right \|$ between the points of the first and second surfaces. We discuss this problem in the Appendix and, using the law of refraction, obtain the inequality of Eq. (27) determining the “permitted” deviation angles. In what follows, we suppose that the functions $\Psi (\mathbf {y})$ and $L(\mathbf {y})$ are defined “correctly” in the sense that the lenses of Eq. (6) focusing to the points $\vec {u}_2(\mathbf {y})$ exist and that the deviation angles satisfy Eq. (27).

4. Supporting quadric method

The supporting quadric method (SQM) [15–21] is widely used in the problems of calculating mirrors and refracting surfaces, which generate prescribed irradiance distributions, and consists in the following. The required continuous irradiance distribution is approximated by a discrete distribution “concentrated” in a finite number of points. The optical surface is represented as a piecewise-smooth surface consisting of fragments of quadrics. Further, the parameters of the quadrics are calculated iteratively from the condition of obtaining the required values of the light flux at the points of the discrete distribution.

In this section, we will consider a version of the SQM for solving the problem of generating prescribed irradiance distribution and wavefront [Fig. 1(a)]. Generally speaking, the rigorous formulation of the SQM requires the introduction of the concept of the so-called weak solution, which allows one to correctly determine the solution of the problem in the case when the ray mapping has discontinuities on a set of measure zero [12,17]. 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}),\, \mathbf {y}\in D$ be approximated by a discrete distribution. To do this, we introduce a mesh with square cells $D_i,\, \, i=\overline {1,N}$ covering the region $D$ and replace $L(\mathbf {y})$ with its following discrete approximation:

The corresponding discrete eikonal function is introduced as $\Psi _i =\Psi (\mathbf {y}_i),\, \, i=\overline {1,N}$. In this case, the vector-function $\vec {p}(\mathbf {y})$ defining the ray directions and the function $l(\mathbf {y})$ defining the second surface are also replaced by discrete sets of values $\vec {p}_i =\vec {p}(\mathbf {y}_i),\, \, i=\overline {1,N}$ and $l_i =l(\mathbf {y}_i)$, $i=\overline {1,N}$. The values $l_i$ will be varied in the method in order to achieve the required flux values $L_i$ at the points $\vec {u}_{2,i} =\vec {u}_2 (\mathbf {y}_i)$, $i=\overline {1,N}$ of the second surface. In this approach, it is assumed that the ray $\vec {p}_i$ will “transfer” the light flux $L_i$ from the second surface to the considered point $\mathbf {y}_i$ of the output plane.

With the discrete approximations introduced above, the surface $R_1$ will be a piecewise-smooth surface consisting of fragments of refractive lenses (6) (ellipsoids) focusing the incident plane beam to the points $\vec {u}_{2,i} =\vec {u}_2 (\mathbf {y}_i)$, $i=\overline {1,N}$. The cross-section of this surface is schematically shown in Fig. 1(b) with highlighted (thick) segments of colored lines showing the ellipsoids focusing to separate points of the second surface. The dashed black line shows the envelope of this family of ellipsoids, which appears in the limiting case when the number of points $N$ tends to infinity.

According to Eq. (8), the function $u_1 (\mathbf {x})$ defining the surface $R_1$ reads as

Now, for fixed values $l_i ,\, \, i=\overline {1,N}$, one can calculate the discrete energy distribution generated at the points $\vec {u}_{2,i}$ with the function $u_1 (\mathbf {x})$ defined by Eq. (12). According to Eqs. (9) and (12), to the point $\vec {u}_2 (\mathbf {y}_i)$, the light flux from the region $C\left (\mathbf {y}_i; l_1, \ldots ,l_N \right )\subset G$ is focused, which is defined by the following condition:

Accordingly, the light flux focused to the point $\vec {u}_{2,i}$ can be calculated as

In practical calculations, the ray tracing approach schematically shown in Fig. 1(b) is commonly used for finding the $L_{\textrm {est},i}$ values (i. e., for calculating the integral in the right-hand side of Eq. (14)). In this case, the incident beam is approximated by a set of $N_\textrm {tr}$ rays going out of the nodes of a certain mesh $\mathbf {x}_j ,\, j=\overline {1,N_\textrm {tr} }$ in the plane $z=0$. These rays “carry” the light flux values $I_j ,\, j=\overline {1,N_\textrm {tr} }$, which are proportional to the irradiance of the incident beam at the mesh nodes. According to Eq. (12), for each ray going out of the point $\mathbf {x}_j$, the index of the lens $i(j)=\mathop {\arg \min }\nolimits _i \left (u_\textrm {lens} (\mathbf {x}_j ;\mathbf {y}_i ,l_i)\right )$ is found, which the ray “meets” first [Fig. 1(b)]. This ray is directed to the corresponding point $\vec {u}_{2,i}$ of the second surface. Then, at the points of the second surface, the light flux values are calculated, which are equal to the sums of the light flux values of the rays arriving to these points:

For the sake of illustration, the calculated light flux values are shown in Fig. 1(b) with integers equal to the number of rays arriving to the points of the second surface.

In the Appendix, it is shown that the function $z=u_\textrm {lens} (\mathbf {x};\mathbf {y},l(\mathbf {y}))$ is monotonic with respect to the parameter $l=l(\mathbf {y})$, and $\partial u_\textrm {lens} / \partial l <0$. This means that with an increase in the $l_i$ value (in the case when the rest values $l_j$, $j\ne i$ are unchanged), the surface $u_\textrm {lens} (\mathbf {x}_j ;\mathbf {y}_i ,l_i)$ “moves down”. In this case, the region $C\left (\mathbf {y}_i ;l_1, \ldots , l_N \right )$ becomes larger and, consequently, the light flux value $L_{\textrm {est},i} =L_{\textrm {est},i} \left (l_1, \ldots , l_N \right )$ increases. Thus, if at a certain point $\vec {u}_{2,i}$ the required light flux $L_{d} (\mathbf {y}_i)=L_i >L_{\textrm {est},i}$, it is necessary to increase the $l_i$ value, and to decrease it in the opposite case. It is straightforward to choose the amount of the increase (or decrease) of $l_i$ proportionally to the difference $\left (L_i -L_{\textrm {est},i} \right )$ between the required and the current calculated values of the light flux [21,22]. Therefore, the $l_i$ values can be calculated iteratively in the following way. Let us denote by $l_{i,n} ,\, i=\overline {1,N}$ the values of $l_i$, obtained at the $n$-th step of the iterative process. Then, at the next [$(n+1)$-th] step, the generated distribution $L_{\textrm {est},i} \left (l_{1,n}, \ldots , l_{N,n} \right ),\, i=\overline {1,N}$ is calculated, and the values $l_{i,n}$ are corrected according to the rule

It is also interesting to note that the SQM of Eq. (16) can be regarded as a method of simple iteration for solving a system of nonlinear equations. Indeed, the problem of calculating the values $l_i ,\, i=\overline {1,N}$, providing the required light flux values $L_i ,\, i=\overline {1,N}$ can be considered as the problem of solving the following system of nonlinear equations:

The system (17) can be rewritten in the following equivalent form:

It should be noted that the simple iteration method is less efficient as compared to the Newton’s method or quasi-Newton methods, which are usually used in design approaches [3–9] based on the solution of an NDE of elliptic type. At the same time, as the computational examples presented in the next section demonstrate, the SQM in the form of Eq. (16) makes it possible to effectively solve the problem of generating light beams with prescribed irradiance distributions and eikonal functions, including an irradiance distribution specified in a doubly-connected region in the form of a square frame (examples in Subsections 5.3 and 5.4). It is important to note that when the irradiance distribution of the output beam is specified in a non-connected region, the methods of [3–9] based on solving an NDE become inapplicable [11,22].

5. Design examples

To illustrate the performance of the method, we designed a number of optical elements generating light beams with prescribed irradiance distributions and wavefronts. In examples 1–3 (Subsections 5.1–5.3), we consider the generation of light beams with spherical wavefronts (converging and diverging). In our opinion, such beams are important in practical problems. In the last example (Subsection 5.4), to illustrate the broad capabilities of the method, we consider the calculation of an optical element that generates a beam with a complex wavefront, providing focusing of the output “square-frame” beam (such a beam has uniform irradiance in a doubly connected region in the form of a square frame) into a circle.

5.1 Generation of a square beam with a uniform irradiance distribution and a converging spherical wavefront

As a first example, we calculated an optical element transforming an incident plane beam with circular cross-section and uniform irradiance

The calculation of the element was carried out for the following parameters: radius of the incident beam $R=1\textrm {~cm}$; the size of the generated square beam $w=2\textrm {~cm}$; output plane $z=f=10\textrm {~cm}$; $z$-coordinate of the focus of the spherical wavefront of the output beam $F=14\textrm {~cm}$, the eikonal value at the point $\mathbf {y}=0$ of the output plane $\Psi _0 =13\textrm {~cm}$; refractive indices of the element and the surrounding medium $n_1 =1.5$ and $n_0 =1$. Let us note that at the chosen parameters $\Psi _0$, $n_0$, $n_1$, the resulting thickness of the optical element at the point $\mathbf {x}=0$ equals $t_0 = 6\textrm {~cm}$.

The SQM described in the previous section allows for the so-called multiscale implementation. This means that the calculation of an element can be carried out in several stages on sequentially refined meshes (i. e., with a sequential increase in the number of points $\mathbf {y}_i ,\, \, i=\overline {1,N}$ constituting the discrete approximations of the required irradiance distribution and wavefront). In this case, the values obtained by interpolating the $l_i$ values calculated on the coarser mesh at the previous stage of the algorithm can be used as an initial approximation for the $l_i$ values defining the second surface of the element on a finer mesh at the next stage. As the calculation results demonstrate, such a multiscale approach provides fast convergence of the iterative process in Eq. (16).

In the present example, the calculation of the surfaces of the optical element was carried out in three stages. At the first stage, the prescribed irradiance distribution (20) and wavefront (21) were approximated by discrete distributions defined on a square mesh containing only $N=N_1 = 20^{2}$ points. At such a small initial number of points, the method of Eq. (16) converges fast (in a few hundred iterations). At the second stage, the prescribed irradiance distribution and wavefront were approximated by discrete distributions on a refined mesh containing $N=N_2 =100^{2}$ points. As the initial values of $l_i$, $i=\overline {1,N_2 }$ on the refined mesh, spline interpolation of the values $l_i$, $i=\overline {1,N_1 }$ calculated at the first stage on a coarse mesh was used. With the refined mesh, 15000 iterations were made. Finally, at the third stage, the finest mesh consisting of $N=N_3 =200^{2}$ points was used. As previously, spline interpolation of the values calculated at the previous stage was used to obtain the initial values of $l_i$, $i=\overline {1,N_3 }$. As in the case of stage 2, 15000 iterations were performed at stage 3. The total computation time of the optical element on a desktop PC (Intel Core i9-7940X CPU) was of about 15 minutes.

In accordance with the SQM description presented in the previous section, at each iteration of each stage, the light flux values at the points $\vec {u}_{2,i} =(\mathbf {y}_i ,f)-\vec {p}_i l_i, i=\overline {1,N}$ were calculated using the ray tracing technique, and then the values $l_i$ were updated using Eq. (16). In the ray tracing method, the incident beam was approximated by $N_\textrm {tr} =4N$ rays defined on a square mesh in the plane $z=0$. Let us also note that to achieve better convergence, at each iteration we decreased the value of the step size $\Delta _n$ by $0.03\%$.

After the calculations, the second surface described by a set of values $l_i$, $i=\overline {1,N}$ was approximated using the method of least squares by a system of two-dimensional B-splines [11,22]. Then, from the obtained second surface, the first surface was calculated using Eq. (8). The calculated surfaces of the optical element are shown in Fig. 2(a). In accordance with the chosen value $\Psi _0 =13\textrm {~cm}$, the thickness of the optical element (distance between the surfaces at the point $\mathbf {x}=0$) amounts to $t_0 =6\textrm {~cm}$. For the designed element, the distance from the second surface to the output plane equals $l(0)=0.53\textrm {~cm}$. This distance is due to the fact that when calculating the element, the initial $l_i$ values at the first iteration were set to be constant: $l_0 =0.5\textrm {~cm}$.

 

Fig. 2. (a) Surfaces of the optical element transforming a circular beam into a square beam with a converging spherical wavefront. (b)–(d) Irradiance distributions generated by the element in different planes. Red and blue lines show the cross-sections of the irradiance distributions.

Download Full Size | PDF

Figures 2(b)–2(d) show the calculated normalized irradiance distributions generated by the designed optical element in the output plane $z=f=10\textrm {~cm}$, in the plane $z=11\textrm {~cm}$, and in the plane $z=12\textrm {~cm}$. In this and the following examples, the generated irradiance distributions were calculated by tracing $N_\textrm {ray} = 10^{7}$ rays. The presented simulation results demonstrate the formation of a square beam. The normalized root-mean-square deviations (NRMSD) of the calculated irradiance distributions from the uniform ones amount to 6.5% ($z=10\textrm {~cm}$), 7.4% ($z=11\textrm {~cm}$), and 8.2% ($z=12\textrm {~cm}$). A linear decrease in the beam size, which occurs with an increase in the distance from the output plane, shows that the wavefront of the generated beam is a converging spherical one. The calculated RMSD of the optical path length in the output plane from the required eikonal function [Eq. (21)] amounts to 14 nm.

5.2 Generation of a square beam with a uniform irradiance distribution and a diverging spherical wavefront

As the next example, we designed an optical element transforming an incident plane beam [Eq. (19)] into a square beam with uniform irradiance [Eq. (20)] and a diverging spherical wavefront. In this case, the focus of the spherical wavefront of the output beam is located in front of the output plane at a certain point $\vec {x}=\left (0,0,F\right )$, where $F<f$. For a diverging spherical wavefront, the eikonal function in the output plane $z=f$ has the form

The element was calculated at $F=6\textrm {~cm}$ and $\Psi _0 =11\textrm {~cm}$. All the rest geometrical and material parameters coincide with the parameters of the previous example. The calculation of the optical element was carried out in three stages and was completely similar to the calculation of the element in example 1. In the present example, the initial values of $l_i$ at the first iteration were chosen to be equal to the constant $l_0 =1\textrm {~cm}$.

The obtained surfaces of the optical element are shown in Fig. 3(a). In accordance with the value $\Psi _0 = 11\textrm {~cm}$ specified in the element design [see Eq. (22)], the thickness of the optical element (distance between the surfaces) equals $t_0 =2\textrm {~cm}$. The distance from the second surface to the output plane equals $l(0)=1.23\textrm {~cm}$. Figures 3(b)–3(d) show the calculated normalized irradiance distributions generated by the element in the output plane $z=f=10\textrm {~cm}$, in the plane $z=12\textrm {~cm}$, and in the plane $z=14\textrm {~cm}$. The presented simulation results demonstrate the generation of a square beam. A linear increase in the beam size with an increase in the $z$ coordinate confirms that the wavefront of the output beam is indeed diverging spherical. The calculated RMSD of the optical path length of the rays in the output plane from the prescribed eikonal function (23) amounts to 11 nm. The NRMSD values of the generated irradiance distributions from constant values equal 6.3%, 4.4%, and 3.4% in the planes $z=10\textrm {~cm}$, $z=12\textrm {~cm}$, and $z=14\textrm {~cm}$, respectively.

 

Fig. 3. (a) Surfaces of the optical element transforming a circular beam into a square beam with a diverging spherical wavefront. (b)–(d) Irradiance distributions generated by the element in different planes. Red and blue lines show the cross-sections of the irradiance distributions.

Download Full Size | PDF

5.3 Generation of a square-frame beam with a diverging spherical wavefront

As a more complex example, we considered the design of an optical element transforming a plane incident beam [Eq. (19)] into a square-frame beam with a diverging wavefront. This beam has a uniform irradiance in a region in the form of a square frame:

The element was calculated for an output beam with the size $w_1 =2\textrm {~cm}$ and $w_2 =1\textrm {~cm}$. The other parameters are the same as in the previous example. The surfaces of the optical element were calculated in several stages with sequential refinement of the meshes similarly to examples 1 and 2. Due to fact that the required distribution has a more complex form and is specified in a doubly-connected region, a mesh of $N=300^{2}$ points was used at the last calculation stage. The resulting total computation time was of about 1 hour.

The calculated surfaces of the optical element are shown in Fig. 4(a). The maximum distance between the surfaces equals $2.38\textrm {~cm}$, and the distance from the second surface to the output plane is $1.14\textrm {~cm}$. Note that the central “square hole” in the second surface, caused by the square-frame shape of the output beam, can be “closed” by an arbitrary surface without affecting the performance of the element. This is due to the fact that the rays refracted by the first surface do not pass through the hole in the second surface. Let us also note that in the vicinity of the origin of coordinates, the first surface has sharp bends 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 vicinity of the point $\mathbf {x}=0$) into the inner contour of the square frame. Such a mapping is discontinuous. The possibility of using the SQM for calculating piecewise-smooth surfaces implementing discontinuous ray mappings is an advantage of the method. In this case, the known methods based on numerical solution of elliptic NDEs [3–9] become inapplicable or turn out to be numerically unstable. In particular, for the application of the “NDE-based” methods in these cases, it is necessary to define the boundary conditions specifying the mapping between the boundary of the incident beam (circle) and the boundary of the target region (square frame). This problem is non-trivial since the circle is simply-connected (1-connected), whereas the square frame is not. Moreover, the boundary of the square frame is not smooth, and, therefore, the construction of a differentiable mapping is impossible.

 

Fig. 4. (a) Surfaces of the optical element transforming a circular beam into a square-frame beam with a diverging spherical wavefront. (b)–(d) Irradiance distributions generated by the element in different planes. Red and blue lines show the cross-sections of the irradiance distributions.

Download Full Size | PDF

Figures 4(b)–4(d) show the calculated normalized irradiance distributions generated by the optical element in the output plane $z=f=10\textrm {~cm}$, in the plane $z=12\textrm {~cm}$, and in the plane $z=14\textrm {~cm}$. The presented simulation results demonstrate the formation of a square-frame beam. The NRMSD values of the calculated irradiance distributions from constant values inside the square frame amount to 8.9% ($z=10\textrm {~cm}$), 8.7% ($z=12\textrm {~cm}$), and 8.1% ($z=14\textrm {~cm}$). A linear increase in the size of the beam with increasing $z$ coordinate confirms that the wavefront of the generated beam is diverging spherical. The calculated RMSD of the optical path length of the rays in the output plane from the prescribed eikonal function (23) equals 17 nm.

5.4 Generation of a square-frame beam with a complex wavefront

In order to demonstrate the capabilities of the SQM in the problem of generating light beams with complex wavefronts, we calculated an optical element transforming a plane beam [Eq. (19)] into a square-frame beam [Eq. (24)] with an eikonal distribution $\Psi (\mathbf {y})$ in the output plane $z=f$ providing focusing of the output beam to a circle in the plane $z=F,\, F>f$.

Let, as in the previous example, the dimensions of the output square-frame beam be equal to $w_1 =2\textrm {~cm}$ and $w_2 =1\textrm {~cm}$. Figure 5(a) shows the eikonal function in the output plane $z=f=10\textrm {~cm}$ calculated from the condition of focusing of the square-frame beam [Eq. (24)] to a circle with the radius $r=1.5\textrm {~cm}$ in the plane $z=F=15\textrm {~cm}$. For the calculation of the eikonal function, a version of the SQM presented by some of the present authors in [21,23] was used. The SQM version of [21,23] is aimed at the problem of calculating an eikonal function of the light field defined in a certain plane from the condition of generating a prescribed irradiance distribution in another plane. As quadrics, eikonal functions of converging spherical beams with the foci at the points of the prescribed discrete distribution (in the present case, at the points of a circle in the plane $z=15\textrm {~cm}$) are used in the SQM of [21,23]. Figure 5(b) shows the irradiance distribution generated in the plane $z=15\textrm {~cm}$ in the case of the eikonal function of Fig. 5(a). This distribution was calculated using the ray tracing technique and corresponds with high accuracy to the required uniform distribution on a circle.

 

Fig. 5. (a) Required eikonal function in the output plane $z=f=10\textrm {~cm}$ calculated from the condition of focusing of the square-frame beam into a circle located in the plane $z=F=15\textrm {~cm}$. (b) Irradiance distribution generated by a square-frame beam with the eikonal function of (a) upon its propagation from the output plane to the plane $z=15\textrm {~cm}$.

Download Full Size | PDF

The calculation of an optical element generating a square-frame beam with the eikonal function of Fig. 5(a) was carried out similarly to examples 1–3. The only difference is that the function used as $\Psi (\mathbf {y})$ is not the eikonal of a spherical beam (converging or diverging), but the function shown in Fig. 5(a). As in the case of examples 1–3, the element was calculated in several stages with a sequential refinement of the computation meshes. At the last stage, a mesh consisting of $N=500^{2}$ points was used. The total computation time for all stages amounted to about 2.5 hours.

The calculated surfaces of the optical element are shown in Fig. 6(a). The maximum distance between the surfaces equals 2.35 cm, and the distance from the second surface to the output plane is 0.81 cm. As in example 3, the first surface has sharp bends in the vicinity of the origin of coordinates and, in this vicinity, resembles a tetrahedral pyramid. Let us remind that the presence of these sharp bends is due to the discontinuous nature of the mapping implemented by the first surface, in which the central part of the surface maps the central part of the incident beam (a vicinity of the point $\mathbf {x}=0$) to the inner contour of the square frame.

 

Fig. 6. (a) Surfaces of the optical element transforming a circular beam into a square-frame beam with an eikonal distribution of Fig. 5(a) providing the subsequent focusing of the output beam into a circle. (b)–(d) Irradiance distributions generated by the element in different planes. Red and blue lines show the cross-sections of the irradiance distributions.

Download Full Size | PDF

Figures 6(b)–6(d) show the calculated normalized irradiance distributions generated by the designed optical element in different planes. The presented results demonstrate the generation of a square-frame beam in the output plane $z=f=10\textrm {~cm}$ [Fig. 6(b)] and the generation of the prescribed irradiance distribution in the form of a circle in the plane $z=15\textrm {~cm}$ [Fig. 6(d)]. The irradiance distribution in an intermediate plane $z=12.5\textrm {~cm}$ [Fig. 6(c)] shows the change in the shape of the beam occurring during its focusing from a square frame into a circle. In the output plane, the NRMSD of the calculated irradiance distribution from the uniform one defined on the square-frame region equals 8.4%. The NRMSD of the calculated irradiance distribution from a uniform distribution defined on a circle in the plane $z=15\textrm {~cm}$ does not exceed 1%. The calculated RMSD of the optical path length of the rays in the output plane from the required eikonal function shown in Fig. 5(a) amounts to 9 nm.

6. Conclusion

We proposed a version of the supporting quadric method for calculating refractive optical elements with two working surfaces that transform an incident beam with a plane wavefront into an output beam with prescribed irradiance distribution and wavefront. Using this method, we calculated optical elements, which transform a circular beam with uniform irradiance distribution into square beams with spherical wavefronts (converging or diverging), as well as into a square-frame beam with a diverging spherical wavefront specified in a doubly-connected region. As the most complex example illustrating the capabilities of the proposed method, an optical element was calculated that generates a square-frame beam with a complex wavefront providing the subsequent focusing of the output beam into a circle.

The presented examples demonstrate high performance of the method. Examples of the generation of square-frame beams show the possibility of using the method for calculating optical elements with continuous piecewise-smooth surfaces.

In this paper, we designed optical elements generating light beams with uniform irradiance distributions in the case of a uniform incident beam. However, the SQM presented above was formulated for a quite general case and, without any modifications, can also be applied in the case of non-uniform distributions of the incident and output beams [22]. The investigation of the SQM performance in the case of complex non-uniform (in particular, raster) irradiance distributions will be a topic of further research.

Due to the reversibility of the light rays, the proposed method can also be used to calculate optical elements transforming an incident beam with an arbitrary wavefront into an output beam with a required irradiance distribution and a plane wavefront. In particular, the practically important problem of forming a collimated optical beam with a prescribed irradiance distribution from a point light source can be solved.

In the general case, when the wavefronts of both the incident and output beams are non-planar, a generalization of the proposed method can be formulated. To obtain this generalized method, one has to rewrite Eqs. (5) and (6) specifying the quadrics. These quadrics should define the lenses that focus the incident beam with a non-planar wavefront to the points of the second surface. For example, with a spherical incident beam, Cartesian ovals will be the quadrics instead of ellipsoids of Eq. (6).