#### Abstract

Beam of light shaping process can be considered ultimate, if both irradiance and wavefront spatial distributions are under control and both can be shaped arbitrarily. In order to keep these two quantities determined simultaneously, it is required to apply at least two powered refractive or reflective surfaces. In this paper, a fully geometric design method of double-freeform beam shapers is discussed briefly. The presented algorithm is based on two stages. First, integrable input-output ray mapping is calculated by the application of the novel GATMA (Geometric Approach to Monge–Ampere equation) method. It allows us to determine the shape of the first freeform surface. Then, according to the condition of constant optical path length between input and output plane, corrected by wavefront phases at those planes, the second surface is determined. GATMA algorithm combines advantages of Monge–Ampere (MA) equation and ray-tracing efficient apparatus. Compared to the state-of-the-art freeform design methods, GATMA does not need to solve MA equation directly but uses this equation as an error function. Such approach makes the computation algorithm simpler and more robust and convergent. The application of the proposed method in a challenging design example of a beam shaper, transforming uniform collimated beam into a beam having a triangular cross section and flat wavefront, is presented as a case study.

#### 1. Introduction

Beam shaping is the process where well-defined input beam is transformed into the desired output beam. It has numerous applications in all the fields where spatial light distribution plays vital role in illumination design [1], solar concentration optics [2], material processing [3], laser active media pumping [4], laser output wavefront correction, semiconductor lasers beam forming [5], lidar [6], military optronics [7], and many others. Beam shaping process can be realized by a variety of optical methods, based on refractive, reflective, diffractive, and adaptive components [8]. A general beam of light is defined by two quantities: irradiance spatial distribution and wavefront (or phase spatial distribution). Beam shaping has to control both of them, and for this reason, beam shapers include at least two essential optical surfaces: one for irradiance and the other for wavefront (or phase) tuning. In some applications, where the only requirement is associated with the desired irradiance that has to be obtained at a fixed distance, a single optical surface can be used. Concerning the state of the art, most of the beam shapers available in the market deal with Gaussian-to-top-hat beam transformation. It is well described, analytically solvable problem addressed by many researchers. Due to rotational symmetry of both Gaussian and top-hat beam, the Gaussian-to-top-hat beam shapers (also called pi-shapers) are composed of traditional (rotationally symmetrical) aspherical components. Manufacturing technology of such components is challenging but available in specialized optical centres. Until recently, the symmetry (rotational or translational) of optical components was the key manufacturability requirement due to the limitations of nonsymmetrical (freeform) optical surfaces fabrication. This has changed however in recent years, which can be easily noticed through the common availability of progressive lenses [9, 10]. Capabilities of precise multiaxis machine tools increased to the level where completely asymmetrical optics can be manufactured [11, 12]. It created new prospects for light shaping possibilities.

Nevertheless, apart from manufacturing, freeform optics is also challenging in terms of design. There are several methods to tackle this problem. Oliker considered freeform surface as the superposition of portions of ellipsoids (or hyperbolics), each of them responsible for one pair of input-output points [13]. The method is conceptually simple; however, its application is reasonable only in case of noncomplicated grids transformation. Another approach was presented by Muschawek [14] who dealt with irradiance shaping from the perspective of wavefronts tailoring. Unfortunately, insufficient computational details were given to reproduce this methodology by other researchers. Another approach is a well-known SMS (simultaneous multiple surfaces) method which was proposed by Minano [15, 16]. The algorithm allows design of single lens with freeform surfaces on both sides. Although it is very efficient method in terms of short distance light transfer optimization, the precision of output wavefront tailoring, especially in case of significant difference between input and output irradiance footprint, is questionable from technological point of view due to limited thickness of reasonably fabricated single lens. In [17], it is demonstrated that freeform lens design can lead to Poisson partial differential equation. The validity of such simplification however is limited to subtle irradiance redistribution. More universal methodology is based on optimal mass transportation and corresponding elliptic Monge–Ampere partial differential equation [18–20] which allows design of single freeform surface capable of creating arbitrary irradiance transformation. In most cases, however, the presented solutions deal only with desired irradiance to be obtained at the fixed output plane. There are few recent papers dealing with the controlled transformation of both irradiance and wavefront [21–24]; however the applied mathematical approaches are based also on optimal mass transportation and corresponding ray mapping.

The presented algorithm differs from all the discussed methods, since it is based purely on ray propagation and ray-distribution-to-irradiance transformations. Each irradiance transformation can be considered in the geometrical optics domain as a specific ray mapping. In 1D case, which can be associated with rotationally or translationally symmetric optical systems, for each transformation there is unique corresponding ray mapping. Its derivation, based in 1D integration, is trivial. In 2D case, the problem becomes more complicated. First of all, each irradiance transformation corresponds to infinite number of ray mappings. It is not difficult to calculate some of those mappings, for example, from differential energy conservation rules; however, the resultant mapping will not probably be integrable. It means that the corresponding refractive or reflective optical surface associated with such mapping will have to be discontinuous and thus hard to manufacture. Integrability condition of optical surface ensures its continuity and is fundamental requirement for fabrication possibility.

In this paper, we deal not only with irradiance, but also with wavefront shaping. By controlling the wavefront, it is possible to dictate how the beam has to propagate beyond the plane of prescribed irradiance. For example, imposing a flat wavefront is associated with beam collimation; the irradiance pattern will remain (in geometrical approximation) unchanged upon propagation. The discussed methodology deals with general beam shapers, composed of two plano-freeform lenses (Figure 1). The first freeform lens realizes integrable ray mapping associated with desired irradiance transformation, and the second lens is responsible for wavefront shaping.

In order to design the first freeform surface, novel approach is proposed. Classical elliptic Monge–Ampere partial differential equation is not solved directly (which sometimes is problematic) but is used to evaluate irradiance error function. According to this philosophy, the described algorithm has got the name GATMA (Geometric Approach to Monge–Ampere equation). Finally, the second freeform surface is obtained by imposing the constant phase shift between input and output wavefront.

#### 2. Materials and Methods

The considered problem can be referred in the configurations presented in Figure 2. Four essential surfaces are considered: input (source) surface *A*, first freeform surface *B*, second freeform surface *C*, and output (target) surface *D*. Additionally, two flat surfaces denoted as *B*′ and *C*′ have been introduced, corresponding to the entrance aperture of the first lens and exit aperture of the second lens, respectively.

Formally, four 2D functions (*I*_{A}, , *I*_{D}, and ), five scalar constraints (three distances, lenses’ depths) are given, and two 2D functions (freeform refractive surfaces of first and second lens) are unknown. The algorithm is divided into three main parts. First, according to the input and output irradiance and input and output wavefront, irradiance distribution at planes *B* and *C* are calculated. It is done by classical ray tracing. Alternatively, diffraction propagation models can be applied. In the second step, optimal ray mapping between planes *B* and *C* is calculated, which is accomplished by irradiance-to-ray distribution driven by the proposed GATMA algorithm. In the last stage, the second freeform surface is designed by forcing for each ray the same optical path length (OPL) between the input and the output wavefront.

##### 2.1. Determination of Optical Field at *B* and *C* Planes

Assume that known irradiance *I*_{A}(*x*_{A}, *y*_{A}) and wavefront (*x*_{A}, *y*_{A}) are associated with *A* plane. Also, the desired output at the *D* plane is predefined: *I*_{D}(*x*_{D}, *y*_{D}) and (*x*_{D}, *y*_{D}). The aim of this step is to calculate irradiance and wavefront at *B* and *C* planes. It can be done in several ways, for example, by the application of well-known diffraction integrals and Fourier methods. In this work, we used purely geometric algorithm, based on ray tracing [25, 26], which makes it very robust and time-efficient. First, irradiance at plane *A* and *D* is transformed into ray distribution, and, according to the wavefront shape, rays are propagated to planes *B* and *C*, respectively. Then, the obtained ray distributions at planes *B* and *C* are transformed into irradiances *I*_{B}(*x*_{B}, *y*_{B}) and *I*_{C}(*x*_{C}, *y*_{C}).

Irradiance is a scalar function, which describes optical power spatial distribution on a considered surface. In ray model, irradiance *I*(*x*, *y*) can be considered as ray density distribution. Assuming that each ray is associated with certain small amount of power Δ*P*, the number of rays *N*(*x*, *y*) impinging onto a small surface element Δ*S* at point (*x*, *y*) can be easily calculated for both *A* and *D* surface using the following formula:

The ray distribution inside elements Δ*S*_{A} and Δ*S*_{D} can be considered as random. As a result, we obtain the discrete distribution of points-rays (*x*_{A,ij}, *y*_{A,ij}) on *A* and (*x*_{D,ij}, *y*_{D,ij}) on *D* surface. Knowing wavefront (*x*_{A}, *y*_{A}) and (*x*_{B}, *y*_{B}), local ray directions on both planes can be represented as unit normal vectors to those wavefronts:

Now, at any ray spot (*x*_{A(D),ij}, *y*_{A(D),ij}), we can create a straight line corresponding to the ray path:where *s* is a scalar parameter which means geometrical path length. It is worth noting that, in case of backward propagation from plane *D* towards plane *C*, *s* has to be negative. On the way between planes *A* and *B* (also *D* and *C*), rays go through the flat input of the first lens (plane *B*′) and the flat output of the second lens (plane *C*′), respectively. In order to determine ray paths **L**_{B′,ij} (s) and **L**_{C′,ij} (s) after refraction on *B*′ and *C*′, ray positions (*x*_{B′,ij}, *y*_{B′,ij}) and (*x*_{C′,ij}, *y*_{C′,ij}) on these planes have to be established. It can be done by the calculation of **L**_{A,ij} and **L**_{D,ij} (*x*, *y*)*-*components, under the condition that *z*-component is equal to *z*_{B′} and *z*_{C′}. This leads to the following formulas:

At each point (*x*_{B′,ij}, *y*_{B′,ij}) or (*x*_{C′,ij}, *y*_{C′,ij}), refraction on a flat surface perpendicular to *z*-axis results in new direction of ray which can be calculated from the vector form of the Snell law:where superscript “(*z*)” indicates the *z*-component of a vector and **z** with a hat is a unit vector along *z*-axis. Accordingly, the paths **L**_{B′(C′),ij} will be given by the following formula:

The distributions of ray spots (*x*_{B,ij}, *y*_{B,ij}) and (*x*_{C,ij}, *y*_{C,ij}) can be obtained from the condition of the above ray lines intersection with *B* and *C* planes, respectively; namely, , which gives

These distributions define the corresponding ray spot density distributions {*N*_{B}; *N*_{C}}, which are obtained by counting how many ray spots are allocated in equal and equally distributed surface cells. Now, the irradiance at both *B* and *C* plane can be obtained by transforming (1) to the following formula:

##### 2.2. Determination of Optimal Ray-Mapping between *B* and *C* Planes

Irradiance transformation on the way between two planes can be equivalently understood in ray-optics sense as ray distribution change or *ray mapping*, which, for each input ray position on the first plane, determines the uniquely corresponding point on the second plane, where this ray should go. Ray mapping can be denoted as a pair of two unknown functions {*x*_{C} = *X*_{C}(*x*_{B}, *y*_{B}), *y*_{C} = *Y*_{C}(*x*_{B}, *y*_{B})}, where (*x*_{B}, *y*_{B}) and (*x*_{C}, *y*_{C}) are ray coordinates at planes *B* and *C*. In the previous section, (*x*_{B,ij}, *y*_{B,ij}) and (*x*_{C,ij}, *y*_{C,ij}) distributions were calculated; however, this was done for irradiance determination at *B* and *C* plane, respectively. Those distributions do not fulfil the integrability condition, which leads to unrealizable optical surface. For this reason, in this section, we briefly describe how to obtain integrable ray mapping.

First of all, irradiance transformation and associated ray mapping have to agree with energy conservation rules:which can also be stated aswhere *J* is the Jacobian 2 × 2 matrix of ray mapping functions. The above equation is mathematically underdetermined, providing numerous possible solutions, which means that theoretically many mappings fulfilling energy conservation rule exist. In practice, however, one has to realize ray mapping by the application of physically realizable optical component associated with smooth (at least of **C**_{1}-class) surface. To meet this requirement, another equation, which is called the integrability condition, has to be met: , where **N** is surface normal vector. Ray mapping associated with integrable optical surface is called integrable ray mapping, or equivalently optimal ray mapping. Among numerous ray mappings fulfilling (10), there is only one integrable ray mapping. At this point, it should be mentioned that the discussed optical problem has a direct connection with OMT (optimal mass transportation) theory; paraxial integrable ray mapping in analogy with optimal mass transport scheme minimizes Monge–Kantorovich integral with quadratic cost function. Following this theory, optimal ray mapping should be a gradient of certain unknown scalar potential *ψ*; namely,

Considering (11) in (10), one obtains the fully nonlinear elliptic partial differential equation of Monge–Ampere type (MA PDE):

It is not analytically solvable; however, there are several applicable numerical methods. This equation was discussed in optical literature, because it allows designing a single freeform surface for arbitrary irradiance shaping on a target plane. In most cases, it is solved directly by finite discretization and the application of Newton method [18, 19]. Alternatively, according to Sulmann [27], it is first transformed into logarithmic parabolic Monge–Ampere equation, and then steady-state solution is searched. In this work, following the purely ray-tracing methodology, the problem described by (12) is not solved mathematically, but geometrically. In the first step, in order to obtain initial ray mapping, we take into account the rays positions at *B* and *C* planes obtained from (7) and denote them as (*x*_{B,ij}, *y*_{B,ij})^{(1)} and (*x*_{C,ij}, *y*_{C,ij})^{(1)}. The corresponding initial ray mapping **M**^{(1)} is thus nonintegrable and does not fulfil (12). This equation can be used however to determine an error function є defined for initial ray positions (*x*_{B,ij}, *y*_{B,ij})^{(1)} and (*x*_{C,ij}, *y*_{C,ij})^{(1)} as follows:

The scalar potential *ѱ*^{(1)}, which appears in (13) is obtained simply from integrating the mapping **M**^{(1)} as it is. In the next step, the error function є^{(1)} is added to , producing slightly changed irradiance distribution at plane *B*, denoted by . Then, the whole loop is repeated *K*-times until є^{(K)} is minimized below certain predefined threshold, which means that the obtained ray mapping corresponding to transformation is “nearly” integrable. Then, by this mapping integration, the optical surface *z*_{B} can be calculated, which is discussed in the following section. It should be underlined that, in the presented method, light distribution at *B* plane is modified, which seems to be incompatible with the statement that this distribution is fixed as an input for the problem. It is true; however, it is also vital that the a.m. modification can be made acceptably small by releasing computation time and working on smaller and smaller differential increments when applying (13).

##### 2.3. Determination of First Freeform Surface

Having the integrable ray mapping determined, we continue with the collection of points (*x*_{B,ij}, *y*_{B,ij}) = (*x*_{B,ij}, *y*_{B,ij})^{(K)}, *i* = 1, …, *M*, *j* = 1, …, *N*, on *B* plane and the corresponding collection (*x*_{C,ij}, *y*_{C,ij}). The size of the collection *M* × *N* should correspond to resolution requirements. Every specific index pair (*i*, *j*) determines one point on *B* plane and the corresponding point on *C* plane. Both points define the associated ray direction after refraction on surface *B*: **r**_{B,ij} = [*x*_{C,ij} − *x*_{B,ij}, *y*_{C,ij} − *y*_{B,ij}, *z*_{C} − *z*_{B}]. Knowing the direction **r**_{B′} of the ray impinging on surface *B*, one can easily determine local slope of the first freeform surface required for a ray to deviate from **r**_{B′} to **r**_{B}. According to the Snell law in vector notation,where **N**_{ij} is a vector which has to be locally normal to the first freeform surface, *z*_{L1}(*x*_{B}, *y*_{B}) at point (*x*_{B,ij}, *y*_{B,ij}). Having the collection of points and normals {(*x*_{B,ij}, *y*_{B,ij}), **N**_{ij}}, one can calculate the surface shape representation: points *z*_{L1,ij} and corresponding 3D point cloud (*x*_{B,ij}, *y*_{B,ij}, *z*_{L1,ij}). Commonly, it is usually done numerically according to Euler scheme; however, in this paper we propose analytical method, based on classical geometry.

Let us assume that the initial point (*x*_{B,11}, *y*_{B,11}, *z*_{B}) on B plane coincides with the first point of freeform surface (*x*_{B,11}, *y*_{B,11}, *z*_{L1,11}). In other words, we state that *z*_{L1,11} = *z*_{B}. It is done without losing generality, since a similar assumption could be made for any *z*_{L1,ij}. The neighboring point in increasing *i*-index direction (*x*_{B,21}, *y*_{B,21}, *z*_{L1,21}) can be obtained as follows. First, we define an equation of the line (ray) crossing the point (*x*_{B,21}, *y*_{B,21}, *z*_{B}) and having the direction **r**_{B′,21}:where *t* is a scalar parameter corresponding to the distance from (*x*_{B,21}, *y*_{B,21}, *z*_{B}). We can also create analytical representation of the surface (*x*, *y*, *z*) perpendicular to **N**_{11} at point (*x*_{B,11}, *y*_{B,11}, *z*_{L,11}):

At this point, it is possible to calculate specific *t* = *t*_{21}, for which line **l**_{21}(*t*_{21}) crosses the surface defined by (16). Combining (15) and (16) leads to the following result:

Identical formula applies for any (*i*, *j*)-index pair and both *i*-index and *j*-index increase directions:

Having *t* determined, it is easy to obtain the *z*_{L1} in neighboring grid points according to the relation , which specifically gives

The above equations can be used to obtain iteratively the whole representation of the first freeform surface *z*_{L1,ij}, *i* = 1, …, *M*, *j* = 1, …, *N*.

##### 2.4. Determination of Second Freeform Surface

Apart from creating the desired irradiance distribution, refraction on the first freeform surface significantly distorts the wavefront. The second freeform surface is introduced to shape the wavefront according to the desired form. The sag *z*_{L2} of this surface can be calculated from the requirement of constant OPL between input and output wavefront for all the rays. Consider the configuration of the beam shaper as presented in Figure 2. It is located between input (*A*) and output (*D*) planes where the wavefront is determined. The first step of the algorithm (par. 2.1) propagates these wavefronts to planes *B* and *C*, respectively. The requirement of constant OPL between wavefronts at *A* and *D* plane is equivalent to the constant OPL between calculated wavefronts at *B* and *C* planes.

At the beginning, the first ray has to be traced and its intersection point with second freeform surface **P**_{C,11} has to be fixed at *C* plane (which translates into *z*_{L2,11} = *z*_{C}). All other intersection points (**P**_{A,11}, **P**_{B′,11}, **P**_{B,11}, **P**_{C′,11}, and **P**_{D,11}) are known from the previous calculations for all rays. **P**_{C,11} coordinates determine the absolute optical path length corresponding to this initial ray propagation through the system:where *n*_{1}, *n*_{2}, and *n*_{3} denote refractive indices of surrounding medium, glass of lenses, and medium inside beam shaper, respectively. To provide wavefront desired transformation, optical path length OPL_{ij} of any other ray (associated with *ij* grid point), traveling between *A* and *D* planes, should be the same as and equal to OPL_{11} provided in (20); additionally, the correction term resulting from input-output wavefront phase differences has to be considered, leading to

It can be further developed into the following equation:

In the above formula, *z*-coordinate of **P**_{C,ij} = (*x*_{C,ij}, *y*_{C,ij}, *z*_{L2,ij}) is the only unknown, so it can be numerically determined at desired precision. **P**_{C,ij} collection creates the point cloud representing the second freeform surface *z*_{L2}(*x*_{C}, *y*_{C}). There is no need for iteration in any predefined sequence, since (22) is pointwise. Additionally, no problems with second freeform surface integrability are faced since the transformed wavefront is assumed to be smooth.

#### 3. Results and Discussion

The presented algorithm was applied to design a beam shaper which can transform a collimated uniform beam of light into uniform triangular irradiance distribution which remains unchanged upon propagation in space. Such property is required to impose a flat wavefront of the output beam. The input beam had a square footprint (size: 20 mm × 20 mm) and a flat wavefront (constant phase spatial distribution). The numerical parameters for this study were the following: *z*_{A} = 0 mm, *z*_{B′} = 25 mm, *z*_{B} = 30 mm, *z*_{C} = 135 mm, *z*_{C′} = 160 mm, *z*_{D} = 500 mm, *n*_{1} = *n*_{3} = 1 (air is both the surrounding medium and the medium between lenses), *n*_{2} = 1.51. The described methodology was implemented in Matlab to calculate the numerical representations of both freeform surfaces (Figure 3).

**(a)**

**(b)**

**(c)**

**(d)**

Additionally, wavefront error at plane *D* was computed, by tracing rays through designed beam shaper. Numerically, wavefront error (represented here by OPD, optical path difference), was calculated by determining, for each ray, the following quantity:

It appeared that wavefront flatness was obtained at very high precision; OPD did not exceed 10^{−12} mm (Figure 4).

In order to verify the designed beam shaper performance, the data was transferred as *grid sag* files into OpticStudio (Zemax). The surfaces were processed as 600 × 600 discrete arrays, which were implemented to create plano-freeform lenses of the beam shaper (Figure 5). Thousands of rays were propagated through the system to simulate its real performance.

The resultant spot diagram was investigated primarily in the output plane (*D* plane) as shown in Figure 6. It can be seen that the required triangular irradiance footprint was successfully obtained.

In order to confirm the wavefront flatness and the corresponding stability of irradiance footprint upon propagation, ray distributions in more distant locations were also analyzed (Figure 7).

**(a)**

**(b)**

In geometrical optics approximation, flat wavefront corresponds to parallel rays. It results in the unchanged beam footprint, no matter how far it propagates. In the case discussed, the distribution shows some distortion if long propagation distance is considered. It indicated that the obtained wavefront is not perfectly flat. The main reason is not the methodology itself but the freeform surface discrete representation finite resolution and limited precision of numerical calculations. The process of transferring Matlab data to OpticStudio, which corresponds to transition from point cloud to solid surface representation of optical element, induces some errors; the surface profile between points is synthesized by simple linear approximation. Taking into account the obtained distortion of triangle beam footprint, it is possible to calculate the corresponding effective net wavefront error. It can be accomplished by the application of the methods used for wavefront reconstruction from Shack–Hartmann wavefront sensor measurements. Specifically, for each ray, knowing its desired position at the distance *z*_{r}: , and its real position *x*_{r}, *y*_{r}, the ray direction error ** δ** (vector) can be calculated from the following formula:

Since a ray is always perpendicular to a wavefront Φ, ray direction error determines wavefront slope (wavefront gradient):

Taking into account (24) and (25), it can be stated thatwhere Φ_{x} and Φ_{y} correspond to partial derivatives ∂Φ/∂*x* and ∂Φ/∂*x*, respectively. Determination of the wavefront from wavefront slopes (wavefront gradient) can be done by several methods. In this research, least-squares minimization scheme on a discrete geometry proposed by Fried [28] was applied. In this approach, both wavefront and its gradient are considered on two-dimensional grid, as presented in Figure 8. Originally, Fried considered a collection of microlenses of a wavefront sensor. In this paper, this approach is modified; instead of (*i*, *j*) microlens, the (*i*, *j*) ray is taken into account.

Assuming perfect ray positions, which form a triangular shape defining a border of beam footprint, wavefront at half grid points inside can be approximated aswhere d*x* and d*y* correspond to the distance between neighboring rays’ perfect positions in *x* and *y* directions. The right-hand sides of the above equations can be written in matrix forms **Γ**_{x}Φ and **Γ**_{y}Φ. Unknown values of Φ at the half-grid points form *N *×* N* (*N*–grid size) array, which can be, by column stacking, transformed into *N*^{2} ×* *1 **Φ-**vector. In the same way, *N *×* N *×* *2 array of known wavefront slopes [Φ_{x}; Φ_{y}] can be stacked into *N*^{2}* *×* *2 **g**-array. Then (27) and (28) can be transformed into matrix form:where **Γ** = [**Γ**_{x}, **Γ**_{y}]^{T}. The least-squares wavefront reconstruction method is reported to provide satisfactory results [29]. In this approach, such solution Φ is searched that minimizes the following term:

Following [29], the minimum-norm solution is given bywhere superscript *T* corresponds to matrix transposition and denotes its pseudoinverse. To apply the discussed method, ray positions at *z* = 100 m were considered as *x*_{r}, *y*_{r} and desired ray positions at *D* plane as . The obtained p–v (peak-to-valley) wavefront error was ∼1.2·10^{−7} m, and rms (root mean square) wavefront error was 0.28·10^{−7} m (Figure 9).

In order to verify the impact of sampling resolution on wavefront flatness, it was increased from 600 × 600 to 1000 × 1000. It resulted in noticeable performance improvement (Figure 10); wavefront error dropped to the level of about 8·10^{−8} m (p–v) and 0.19·10^{−8} m (rms), but at the cost of computational time, which increased about 7 times. Concerning long range propagation of optical beams, one should take into consideration not only the quality of beam shaping optics, but also the atmospheric impact on wavefront distortion caused by inevitable turbulences (nonuniform and dynamic fluctuations of refractive index spatial distribution).

Following both Kolmogorov and non-Kolmogorov models, even moderate conditions lead to significant wavefront angle-of-incidence (AOI) fluctuations and higher order wavefront errors [30]. For this reason, it is always reasonable to keep in mind that, in designing a beam shaper for long-range open-space applications, optics is not limited to lenses but also includes a long column of atmosphere which can be represented as a random phase screen. Sometimes, it is an atmospheric factor which sets the limit of net optical performance, and there is no point in any further improvements of lenses.

#### 4. Conclusions

The paper presents a systematic approach of designing a double-freeform-lens beam shaper capable of controlling the outgoing beam completely, by imposing both desired irradiance and wavefront simultaneously. Most reported methods of freeform lens application in the field of beam transformation deal only with irradiance shaping at a single predefined distance. By controlling also a wavefront, it is possible to define the beam behavior in front of and beyond this distance. Only few papers discuss this issue and propose mathematically complicated methods to solve this problem. In this work, a new “fully geometrical” computational scheme is proposed. As such, it is convergent and relatively easy to follow, also for those who are not familiar with solving partial differential equations. In this approach, elliptic nonlinear Monge–Ampere partial differential equation is also present; however, it is not solved but serves as a merit (error) function of the trial ray tracings.

The proposed method itself deals with input and output planes where the required irradiance and wavefront distributions have to be defined. Additionally, optical designer has to decide where to locate the first and the second freeform lens with respect to those planes. Light propagation from the input plane to the first lens and backpropagation from the output plane to the second lens are modeled by the application of ray-tracing and ray-density-irradiance transformations. The essence of the proposed algorithm lies in the novel approach to find integrable ray mapping between both lenses, which allows us to obtain the geometry of the first freeform surface, responsible for irradiance shaping. The second lens is calculated from the condition of constant input-output OPL corrected by wavefronts’ phase relation.

The algorithm performance was verified by its implementation in a challenging problem to design a beam shaper transforming square top-hat beam into a beam having triangular irradiance footprint, which additionally does not change upon propagation (requirement for perfectly flat wavefront). The designed two-lens system was verified independently in Matlab and OpticStudio (Zemax). In the first case, theoretical wavefront error induced by the system was numerically calculated in the form of OPD map. It revealed excellent performance with OPD fluctuations on the order of 10^{−13} m. Then, the design was transferred to OpticStudio in order to perform extensive ray-tracing analysis. It allowed us to prove the algorithm efficiency; the main goal of the case study challenge (square-to-triangle-footprint combined with flat-to-flat-wavefront beam transformation) was successfully achieved. Additionally, beam irradiance distribution was analyzed far beyond the output plane. According to the flat wavefront property, beams should remain unchanged upon propagation (in geometrical optics approximation). In the presented case study, however, minor beam distortion was observed in OpticStudio simulations. The primary reason for this situation was finite resolution of lenses representations and data transfer process from Matlab to OpticStudio. In the presented research, lenses were designed in Matlab and then exported to OpticStudio in the form of discrete representations of continuous surfaces. OpticStudio applies its own methods to reconstruct such surfaces. To analyze the effect of minor wavefront distortion quantitatively, ray positions predicted by OpticStudio in large distance (1000 m) were compared with their perfect locations (determined by uniform triangular footprint). This allowed us to assess the wavefront residual error, which was done by the application of computation method known from Shack–Hartmann wavefront sensor data analysis. It revealed that the residual rms wavefront error was as small as 0.28·10^{−7} m. In order to prove the resolution impact, it was increased from 600 × 600 to 1000 × 1000 and the wavefront error analysis repeated. It resulted in more than one order of magnitude improvement; rms wavefront error dropped to 0.19·10 ^{−8} m.

#### Data Availability

The data used to support the findings of this study are included within the article.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

This research was supported by internal resources of the Institute of Optoelectronics, Military University of Technology.