Supporting Quadric Method for Designing Freeform Mirrors That Generate Prescribed Near-Field Irradiance Distributions

: We consider a version of the supporting quadric method for designing freeform mirrors that generate prescribed irradiance distributions in the near ﬁeld. The method is derived for a general case of an incident beam with an arbitrary wavefront. As an example, for a practically important special case of a plane incident wavefront, we design a freeform mirror that generates a complex-shaped uniform irradiance distribution in the form of the abbreviation “IPSI” on a zero background. The designed mirror is fabricated and qualitatively investigated in a proof-of-concept optical experiment. The experimental results conﬁrm the correctness of the proposed approach and demonstrate the manufacturability of the mirrors designed using the considered method.


Introduction
The problem of designing a reflecting or refracting optical surface that can generate a required irradiance distribution in a certain domain belongs to the inverse problems of nonimaging optics and is extremely complex. In the geometrical optics approximation, this inverse problem can be reduced to finding a solution to a nonlinear differential equation (NDE) of an elliptic type [1][2][3][4][5][6][7][8][9][10]. Although several finite difference methods have been proposed for solving such NDEs [1][2][3][4][5][6][7][8][9][10], the design of optical elements using this approach is sophisticated and has significant limitations. In particular, the formulation of the problem of the calculation of an optical element as an NDE assumes that the optical surfaces of the element are smooth. The requirement of smoothness limits the class of the irradiance distributions that can be generated by the element. For example, an optical element with smooth surfaces cannot generate an irradiance distribution defined in a disconnected region or in a region with non-smooth boundaries [11][12][13][14].
One of the methods widely used for calculating the reflection and refraction of optical surfaces is the supporting quadric method (SQM) [15][16][17][18][19][20][21][22]. The main advantages of the SQM are its versatility, simplicity, and the capability to make calculations relating to optical surfaces that generate "discontinuous" irradiance distributions defined in nonsimply connected regions with complex non-smooth boundaries. With this method, the required irradiance distribution is approximated by a discrete distribution defined on a finite set of points. Then, the optical surface is represented as a set of segments of stigmatic surfaces (referred to as quadrics), focusing the incident beam on the points of the discrete distribution. Depending on the design problem, paraboloids, ellipsoids, hyperboloids, or more complex surfaces, e.g., Cartesian ovals [19], are used as these stigmatic surfaces. The Photonics 2022, 9,118 2 of 12 parameters of the stigmatic surfaces (quadrics) are calculated iteratively from the condition of generating the prescribed discrete distribution. In the general case, the convergence of this iterative method remains a research issue. In some particular cases, when the considered inverse problem of calculating the optical surface can be formulated as a Monge-Kantorovich mass transportation problem (MTP), the supporting quadric method can be regarded as a gradient descent method for minimizing a certain convex function [22]. In this case, it can be stated that the method converges. Such inverse problems include, for example, the problems of designing the refracting surfaces [11] and mirrors [12,23,24] operating in the far field in the case of a spherical or plane incident beam, the problems of collimated beam shaping [22,25], and the problems of calculating the eikonal function of the light field providing the generation of a required irradiance distribution [26]. In the problems, which cannot be reformulated as MTPs, and in particular in the problems of generating a required irradiance distribution in the near field using mirrors or refracting surfaces, the equivalence of the SQM to the gradient descent method has not been proven, and the theoretical justification of its convergence remains an open question. Nevertheless, in solving practical problems, the method usually demonstrates good performance.
It is important to note that the form of the quadrics, the representation of the optical surface, and the iterative algorithm of updating the quadric parameters differ significantly for different problems. In the present work, we consider for the first time a version of the SQM for calculating a reflecting surface that generates a prescribed irradiance distribution in the near field for an incident beam with an arbitrary wavefront, and provide an equation for the stigmatic surface and a detailed description of the method. The good performance of the method is illustrated with a designed example of a mirror that generates a discontinuous irradiance distribution in the form of framed letters "IPSI". The calculation results are confirmed by a proof-of-concept optical experiment.

Problem Statement
Let us consider a three-dimensional space E 3 with the coordinates (x 1 , x 2 , z), in which a light beam with a wavefront W is propagated. The wavefront is defined by a vector function → w(u) = (w 1 (u), w 2 (u), w 3 (u)), where u = (u 1 , u 2 ) ∈ G are certain curvilinear coordinates and G is a domain, in which the beam has non-zero irradiance E(u). Without the loss of generality, let us assume that the beam propagates in a medium with a unit refractive index. The directions of the rays of the light beam coincide with the directions of the unit normal vectors to the wavefront, which are defined by the relation w(u)/∂u i , i = 1, 2, and the symbol "×" denotes the cross-product.
Let the light beam impinge on a mirror with the surface Q ( Figure 1). We define the surface Q through a function l(u) that describes the distance from the wavefront surface to the mirror along the ray propagation directions in the following form: Let the ray going out of the point of the wavefront W with the coordinates u ∈ G be reflected by the mirror and come to a certain point x = (x 1 , x 2 ) of the plane z = 0 ( Figure 1). This defines a ray mapping x = γ(u) that is implemented by the mirror. Note that the directions of the reflected rays are determined by Snell's law in the form q u 2 is a unit normal vector to the mirror surface, and the symbol "·" denotes the dot product. Taking into account Equation (2), the mapping x = γ(u) can be written as where t(u) = −q 3 (u)/r q,3 (u) is the distance from the mirror surface to the plane z = 0 along the direction of the ray reflected by the mirror.
Photonics 2022, 9, x FOR PEER REVIEW 3 of 12 Figure 1. Geometry of the mirror design problem.
Let the ray going out of the point of the wavefront W with the coordinates G  u be reflected by the mirror and come to a certain point Figure 1). This defines a ray mapping ( )   x u that is implemented by the mirror. Note that the directions of the reflected rays are determined by Snell's law in the form is a unit normal vector to the mirror surface, and the symbol "  " denotes the dot product. Taking into account Equation (2), the mapping ( )   x u can be written as  The mapping x = γ(u) defines the irradiance distribution L(x) generated in the plane z = 0. This distribution can be found using the light flux conservation law, where J γ (u) is the Jacobian of the mapping γ. For the further derivations, let us intro- The integral of this function over a domain U ⊂ G defines the light flux through this domain on the surface of the wavefront W as Φ(U) = U I(u)du. Using the introduced function, the light flux conservation law can also be written in the integral form: where B is an arbitrary subset in the plane z = 0. In the present work, we consider an inverse problem relating to the calculation of a mirror surface Q defined through a function l(u) from the condition of generating a required irradiance distribution L(x), x ∈ D in the domain D of the plane z = 0 ( Figure 1).
In the case of a smooth surface Q, Equations (3) and (4) reduce the calculation of the function l(u) to solving a nonlinear partial differential equation of an elliptic type (an equation of Monge-Ampère type) [6][7][8]. In what follows, we consider a method for solving the formulated inverse problem, which does not assume the smoothness of the surface Q and does not require solving the mentioned NDE.

Envelope Representation of the Mirror Surface
The mirror surface can be represented as an envelope of a family of stigmatic surfaces focusing the incident light beam to the points of the target domain D [12,[22][23][24][25]27]. Similar to Equation (1), the equation of a reflecting stigmatic surface focusing the incident beam to the point x = (x 1 , x 2 ) ∈ D of the plane z = 0 can be written as where l st (u; x) is the distance from the wavefront surface to the stigmatic surface along the ray propagation direction. Let us note that the coordinates of the point x = (x 1 , x 2 ) ∈ D in Equation (6) are considered to be parameters. Let us define the function l st (u; x) in Equation (6) from the condition of a constant path length F(x) from the wavefront W of the incident beam to the focus Since the function l st (u; x) also depends on the value F(x) defining the optical path length for a given point x = (x 1 , x 2 ), we will further write it in the form l st (u; x, F(x)). From Equation (7), the function l st (u; x, F(x)) can be obtained as The equation of the envelope surface for the family of surfaces → q st (u; x) with respect to the parameter x ∈ D is defined by Equation (6) and the following two equations [27,28]: , and the parentheses denote the mixed product of the vectors. Substituting Equation (6) into Equation (9) and performing simple transformations, we obtain the envelope equation in the form where l st,x i = ∂l st (u; x, F(x))/∂x i , i = 1, 2. Note that the fulfillment of the conditions Photonics 2022, 9, 118 5 of 12 is sufficient for the fulfillment of the envelope Equation (10). In order to obtain an explicit expression for the envelope, it is necessary to express from Equation (11) the variables (x 1 , x 2 ) as functions of the variables (u 1 , u 2 ) and substitute the obtained expressions to Equation (6). Note that at a fixed value of u = (u 1 , u 2 ) ∈ G, Equation (11) corresponds to a critical point of the function l st (u; x, F(x)) with respect to the variables x 1 , x 2 . As in the works [12,14,22,[29][30][31], in what follows, we will consider an envelope of a special kind, for which Equation (11) defines not only a critical point, but a minimum point. In this case, the functions x(u) = (x 1 (u), x 2 (u)) satisfying Equation (11) will have the form From Equations (12), (11), and (6), it follows that the envelope equation has the form where It is easy to understand [11][12][13][14]22,29,31] that Equation (12) defines the ray mapping x = γ(u) implemented by the envelope surface. Note that the ray mapping and the envelope are completely defined by the function F(x), which determines the optical path length from the wavefront W to the points of the domain D. This makes it possible to formulate the problem of designing a mirror that generates a prescribed irradiance distribution L(x), x ∈ D in the plane z = 0 as a problem of calculating such a function F(x), so that the corresponding ray mapping of Equation (12) satisfies the light flux conservation law defined by Equation (4) or Equation (5).

Supporting Quadric Method
In this subsection, we consider a version of the supporting quadric method (SQM) for the solution of the inverse problem of calculating a mirror that generates a prescribed irradiance distribution ( Figure 1). Strictly speaking, the rigorous formulation of the SQM requires the introduction of the concept of the so-called weak solution, which makes it possible to correctly define the solution of the problem in this case, in which the ray mapping has discontinuities on a set of measure zero [15,17,29,31]. However, for the sake of simplicity, let us focus on the practical side of the method.
In the SQM, the prescribed continuous irradiance distribution L(x), x ∈ D has to be approximated by a discrete distribution. Therefore, let us assume that in the domain D, a rectangular grid containing N cells D i , i = 1, N is introduced, with the centers at the points x i , i = 1, N. As a discrete approximation of the continuous distribution L(x), x ∈ D, let us consider the following discrete distribution: where δ(x) is the Dirac delta function, and L i =  (15), the function F(x) becomes a set of numbers F = {F 1 , . . . , F N }, which we will refer to as the weights of the points In this case, the envelope surface of Equations (13) and (14) turns into a piecewise smooth surface consisting of N segments of the stigmatic surfaces of Equations (6) and (8), focusing the incident beam to the points This surface is defined by Equation (13), where the function l(u) reads as The ray mapping implemented by this surface is defined by the equation Since for a discrete distribution, the surface of Equations (13) and (16) is completely defined by the set of weights F = {F 1 , . . . , F N }, the problem of generating a given discrete distribution of Equation (15) can be considered to be the problem of finding weights Let us explain, how the light flux values at the points x i , i = 1, N are related with the weights F = {F 1 , . . . , F N }. For this, let us introduce the concept of the generalized Voronoi cell to define the "apertures" of the stigmatic surfaces in the domain G. Under generalized Voronoi cells with the weights F = {F 1 , . . . , F N }, we will understand the subdomains (17), it follows that the cells C F (x i ) can be defined by the following condition: a point u ∈ G lies inside the cell C F (x i ) if, and only if, The light flux Φ i focused on the point x i is calculated as an integral over the corresponding generalized Voronoi cell: According to Equation (19), the problem of finding the weights F = {F 1 , . . . , F N } providing the required light flux values L i , i = 1, N at the points x i , i = 1, N can be considered to be the problem of solving the following system of nonlinear equations: Let us rewrite this system as Using the simple iteration method for solving the system (21), we obtain the following iterative algorithm for the calculation of the weight values: where the superscript k denotes the iteration number, and ε > 0 is the step of the iterative method. Let us explain that the algorithm of Equation (22) updates the weights F = {F 1 , . . . , F N } in the "right direction". Assume that all the weights except the value F i are fixed. Let, at the k-th step of the iterative algorithm, the light flux Φ i (F) generated at the point x i be greater than the required one L i . In this case, the algorithm of Equation (22) increases the weight F i . Below, we demonstrate that ∂l st,i /∂F i ≥ 0 and, therefore, for all points u ∈ G the value l st,i (u; x i , F i ) increases. As a result, the inequalities in Equation (18) hold for a smaller set of points u ∈ G, i.e., the cell C F (x i ) decreases and, correspondingly, the generated light flux Φ i (F) also decreases. Similarly, one can show that if the light flux Φ i (F) generated at the point x i is smaller than the required one L i , then the algorithm of Equation (22) will decrease the weight F i . In this case, the size of the cell C F (x i ) and, correspondingly, the generated light flux Φ i (F) will increase.
Let us now prove the following inequality used above: By expanding the right-hand side of Equation (23), we obtain Let us consider the last two terms of the numerator: where ϕ is the angle between the vectors → n w (u) and To conclude this section, let us note that in a practical implementation of the iterative algorithm of Equation (22), the ray-tracing method is used for calculating the values Φ i (F). In this case, the incident beam is approximated by a set of N tr rays going out of the nodes of a certain grid u j , j = 1, N tr on the surface of the wavefront W. These rays "carry" the light flux values I j , j = 1, N tr , which are proportional to the values of the irradiance of the incident beam at the nodes of the grid. According to Equations (16) and (17), for each ray going out of the point u j , the index of the surface i(j) = argmin i (l st,i (u; x i , F)) is found, which the ray "meets" first. This ray is directed to the corresponding point x i of the generated distribution. Then, at the points x i , the light flux values are calculated as the sums of the light fluxes of the rays arriving at these points.

Results and Discussion
The SQM for designing mirrors considered in the previous section was implemented in a practically important case of an incident beam with a plane wavefront. In the example considered below, the wavefront surface W is the plane where → n w = −(1, 0, 1)/ √ 2 and → x 0 = (10 mm, 0, 40 mm). Equation (26) can be written in the parametric form: where → k 1 = (1, 0, −1) · 1/ √ 2, → k 2 = (0, 1, 0) are the unit vectors of the Cartesian coordinate system u = (u 1 , u 2 ) in the plane W with the origin at the point → x 0 (Figure 2). It can be easily understood that in the case of a plane wavefront of the incident beam, the stigmatic surfaces of Equations (6) and (8) that focus the incident light beam to the points of the target domain correspond to paraboloids of revolution. The axes of the paraboloids are parallel to the normal vector → n w to the wavefront, whereas their foci lie in the target domain D. Note that in the plane wavefront case, the only simplification in the general "envelope representation" of the mirror surface (Equations (13) and (14)) consists of the fact that the normal vector → n w (u) does not depend on the point u. , u u  u in the plane W with the origin at the point 0 x  (Figure 2). It can be easily understood that in the case of a plane wavefront of the incident beam, the stigmatic surfaces of Equations (6) and (8) that focus the incident light beam to the points of the target domain correspond to paraboloids of revolution. The axes of the paraboloids are parallel to the normal vector w n  to the wavefront, whereas their foci lie in the target domain D. Note that in the plane wavefront case, the only simplification in the general "envelope representation" of the mirror surface (Equations (13) and (14)) consists of the fact that the normal vector ( ) w n u  does not depend on the point u .
In accordance with the version of the SQM considered in the previous section, the required irradiance distribution ( ), L D  x x (Figure 2) was approximated by a discrete distribution containing 29403 N  points. The points were located at the nodes of a square grid (grid step 0.09   mm) covering the domain D. The calculation of the mirror surface was carried out using the iterative algorithm of Equation (22). At each iteration of the method, the light flux values at the points of the discrete distribution were Figure 2. Geometry of the problem of designing a mirror for generating an irradiance distribution in the form of framed letters "IPSI" for an incident beam with a plane wavefront. The dimensions of the input beam, freeform mirror, and the generated distribution correspond to the considered example up to a scale factor.
Let the incident beam be defined in a rectangular domain G = { u| |u 1 | < 12.5 mm, |u 2 | < 8.5 mm} and have a constant irradiance distribution E(u) = E 0 , u ∈ G (Figure 2). For this incident beam, we designed a mirror generating in the plane z = 0 an irradiance distribution L(x), x ∈ D in the form of framed letters "IPSI" (the abbreviation of the Image Processing Systems Institute) on a zero background (Figure 2). The outer boundary of the frame coincides with the boundary of the square region D = {x| |x 1 − 32 mm| < 20 mm, |x 2 | < 20 mm }.
In accordance with the version of the SQM considered in the previous section, the required irradiance distribution L(x), x ∈ D (Figure 2) was approximated by a discrete distribution containing N = 29, 403 points. The points were located at the nodes of a square grid (grid step ∆ = 0.09 mm) covering the domain D. The calculation of the mirror surface was carried out using the iterative algorithm of Equation (22). At each iteration of the method, the light flux values at the points of the discrete distribution were calculated using the ray-tracing technique, and then the weight values F i were corrected using Equation (22). As the initial weight values for starting the iterative process, constant values were chosen: The F 0 value corresponds to the optical path length of the "central" ray coming out of the point → x 0 reflected from the plane x 1 = 0. In this case, the initial mirror surface is given by Equations (13), (16), and (27).
In the ray-tracing method used at each iteration for calculating the light flux values Φ i , the incident beam was approximated by 3N rays defined on a square grid. In the calculation, about 35,000 iterations were performed. The total computation time on a modern desktop PC (Core i9-7940X CPU) was about three hours.
The three-dimensional image of the designed mirror is shown in Figure 2. On the surface of the mirror, "sharp bends" (derivative discontinuities) can be seen. The presence of the sharp bends is caused by the discontinuous character of the ray mapping, corresponding to the generated "discontinuous" irradiance distribution in the form of separate letters and a frame on a zero background. The derivative discontinuities are additionally shown with blue lines in Figure 3a. This figure shows the calculated mirror as the function x 1 = x 1 (x 2 , z) describing the "depth" of the mirror surface along the x 1 axis. Figure 3b shows the irradiance distribution generated by the mirror in the plane z = 0. This distribution was calculated using the commercially available ray-tracing software TracePro [32], and with a good accuracy coincides with the required distribution. The normalized root mean square deviation of the calculated distribution (Figure 3b) from the required one ( Figure 2) amounts to 8.1%.
Photonics 2022, 9, x FOR PEER REVIEW 9 of 12 calculated using the ray-tracing technique, and then the weight values i F were corrected using Equation (22). As the initial weight values for starting the iterative process, constant values were chosen: . The 0 F value corresponds to the optical path length of the "central" ray coming out of the point 0 x  reflected from the plane 1 0 x  . In this case, the initial mirror surface is given by Equations (13), (16), and (27).
In the ray-tracing method used at each iteration for calculating the light flux values i  , the incident beam was approximated by 3N rays defined on a square grid. In the calculation, about 35,000 iterations were performed. The total computation time on a modern desktop PC (Core i9-7940X CPU) was about three hours. The three-dimensional image of the designed mirror is shown in Figure 2. On the surface of the mirror, "sharp bends" (derivative discontinuities) can be seen. The presence of the sharp bends is caused by the discontinuous character of the ray mapping, corresponding to the generated "discontinuous" irradiance distribution in the form of separate letters and a frame on a zero background. The derivative discontinuities are additionally shown with blue lines in Figure 3a. This figure shows the calculated mirror as the function   1 1 2 , x x x z  describing the "depth" of the mirror surface along the 1 x axis. Figure 3b shows the irradiance distribution generated by the mirror in the plane 0 z  . This distribution was calculated using the commercially available ray-tracing software TracePro [32], and with a good accuracy coincides with the required distribution. The normalized root mean square deviation of the calculated distribution [ (Figure 3b]) from the required one ( Figure 2) amounts to 8.1%.
(a) (b) Figure 3. (a) Surface of the mirror generating an irradiance distribution in the form of framed letters "IPSI" shown as "depth" function along the 1 x axis. Blue lines show the "sharp bends" of the surface (derivative discontinuities). (b) Normalized irradiance distribution generated by the mirror calculated using TracePro.
The presented example demonstrates one of the advantages of the proposed method, which consists of the capability to calculate continuous piecewise smooth surfaces for implementing discontinuous ray mappings. In this case, the known methods based on Figure 3. (a) Surface of the mirror generating an irradiance distribution in the form of framed letters "IPSI" shown as "depth" function along the x 1 axis. Blue lines show the "sharp bends" of the surface (derivative discontinuities). (b) Normalized irradiance distribution generated by the mirror calculated using TracePro.
The presented example demonstrates one of the advantages of the proposed method, which consists of the capability to calculate continuous piecewise smooth surfaces for implementing discontinuous ray mappings. In this case, the known methods based on numerical solution of NDEs of an elliptic type [1][2][3][4][5][6][7][8][9][10] are either inapplicable or numerically unstable.
The designed mirror was fabricated using the milling technique with a three-axis machining center, Haas Minimill. As a blank for the mirror, a silver plate was used. The photograph of the two fabricated mirrors is shown in Figure 4. In the right-hand side of the photograph, a mirror after the milling process is shown, whereas the left-hand side shows a mirror after the polishing process, which was applied in order to remove the cutter marks.

unstable.
The designed mirror was fabricated using the milling technique with a three-axis machining center, Haas Minimill. As a blank for the mirror, a silver plate was used. The photograph of the two fabricated mirrors is shown in Figure 4. In the right-hand side of the photograph, a mirror after the milling process is shown, whereas the left-hand side shows a mirror after the polishing process, which was applied in order to remove the cutter marks.   Figure 4b shows the irradiance distribution generated by the mirror when illuminated by a handheld flashlight. The dimensions and shape of the generated distribution qualitatively correspond to the calculated distribution shown in Figure 3b and confirm the correctness of the performed calculations. The blurring of the letters and artifacts in the generated distribution are caused by the fabrication imperfections and errors in the positioning of the illuminating beam, as well as its deviation from the beam with a plane wavefront and uniform irradiance distribution used in the mirror design. However, the authors believe that the presented example demonstrates the "robustness" of the mirror designed using the proposed version of the SQM to the fabrication errors and the variations in the incident beam. This makes it possible to use such mirrors in "out-of-laboratory conditions", for example, when creating jewelry exhibiting interesting optical effects. In particular, the mirrors shown in Figure 4a are prototypes for cufflinks.

Conclusions
We proposed a version of the supporting quadric method for calculating a mirror that generates a prescribed irradiance distribution in the near field for an incident beam with an arbitrary wavefront. In the method, the mirror surface is represented as a set of segments of stigmatic surfaces focusing the incident beam to the points of the required distribution. The parameters of the stigmatic surfaces are calculated using an iterative method, which can be considered to be a simple iteration method for solving a system of nonlinear equations.
Using the developed method, a mirror generating an irradiance distribution in the form of framed letters "IPSI" on a zero background was designed for a practically important particular case of an incident beam with a plane wavefront. This example demonstrates the possibility of using the method for designing mirrors with continuous piecewise smooth surfaces for implementing discontinuous ray mappings. The numerical simulation results obtained using the commercial software TracePro demonstrate the good quality of the generated distribution. The designed mirror was fabricated using the  Figure 4b shows the irradiance distribution generated by the mirror when illuminated by a handheld flashlight. The dimensions and shape of the generated distribution qualitatively correspond to the calculated distribution shown in Figure 3b and confirm the correctness of the performed calculations. The blurring of the letters and artifacts in the generated distribution are caused by the fabrication imperfections and errors in the positioning of the illuminating beam, as well as its deviation from the beam with a plane wavefront and uniform irradiance distribution used in the mirror design. However, the authors believe that the presented example demonstrates the "robustness" of the mirror designed using the proposed version of the SQM to the fabrication errors and the variations in the incident beam. This makes it possible to use such mirrors in "out-of-laboratory conditions", for example, when creating jewelry exhibiting interesting optical effects. In particular, the mirrors shown in Figure 4a are prototypes for cufflinks.

Conclusions
We proposed a version of the supporting quadric method for calculating a mirror that generates a prescribed irradiance distribution in the near field for an incident beam with an arbitrary wavefront. In the method, the mirror surface is represented as a set of segments of stigmatic surfaces focusing the incident beam to the points of the required distribution. The parameters of the stigmatic surfaces are calculated using an iterative method, which can be considered to be a simple iteration method for solving a system of nonlinear equations.
Using the developed method, a mirror generating an irradiance distribution in the form of framed letters "IPSI" on a zero background was designed for a practically important particular case of an incident beam with a plane wavefront. This example demonstrates the possibility of using the method for designing mirrors with continuous piecewise smooth surfaces for implementing discontinuous ray mappings. The numerical simulation results obtained using the commercial software TracePro demonstrate the good quality of the generated distribution. The designed mirror was fabricated using the milling technique and was qualitatively investigated in a proof-of-concept optical experiment. The experimental results confirm the correctness of the proposed approach and demonstrate the manufacturability of the mirrors designed using the developed method. Funding: This work was supported by the Russian Science Foundation (project 18-19-00326, development of the supporting quadric method in the general case of an incident beam with an arbitrary wavefront, design and fabrication of the mirror using the milling technique) and by the Ministry of Science and Higher Education of the Russian Federation (state assignment to the FSRC "Crystallography and Photonics" RAS (agreement 007-GZ/Ch3363/26), technical support in the CAD design and additional processing (polishing) of the mirror surface; and state assignment for research to Samara University (laboratory "Photonics for Smart Home and Smart City", project FSSS-2021-0016), software implementation of the supporting quadric method).

Data Availability Statement:
The data that support the presented results are available from the corresponding author upon reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest.