# Simplified Geometric Approach to Freeform Beam Shaper Design

Authors: Jacek Wojtanowski; Tadeusz Drozd
Journal: International Journal of Optics
 (2020)

---

## 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.


---

## Body

## 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 (Figure1). The first freeform lens realizes integrable ray mapping associated with desired irradiance transformation, and the second lens is responsible for wavefront shaping.Figure 1
Beam shaper composed of two plano-freeform lenses for arbitrary input/output beams (freeform surfaces: 1 realizes irradiance transformation; 2 adjusts wavefront).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 Figure2. 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.Figure 2
Setup and symbol definitions of the designed double plano-freeform lens beam shaper.Formally, four 2D functions (IA, wA, ID, and wD), 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 atB and C Planes

Assume that known irradianceIA(xA, yA) and wavefront wA(xA, yA) are associated with A plane. Also, the desired output at the D plane is predefined: ID(xD, yD) and wD(xD, yD). 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 IB(xB, yB) and IC(xC, yC).Irradiance is a scalar function, which describes optical power spatial distribution on a considered surface. In ray model, irradianceI(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:(1)NAxA,yA=IAxA,yA⋅ΔSAΔP,NDxD,yD=IDxD,yD⋅ΔSDΔP.The ray distribution inside elements ΔSA and ΔSD can be considered as random. As a result, we obtain the discrete distribution of points-rays (xA,ij, yA,ij) on A and (xD,ij, yD,ij) on D surface. Knowing wavefront wA(xA, yA) and wB(xB, yB), local ray directions on both planes can be represented as unit normal vectors to those wavefronts:(2)rA,ij=−∂wA/∂x,−∂wA/∂y,1∂wA/∂x2+∂wA/∂y2+1x=xA,ij,y=yA,ij,rD,ij=−∂wD/∂x,−∂wD/∂y,1∂wD/∂x2+∂wD/∂y2+1x=xD,ij,y=yD,ij.Now, at any ray spot (xA(D),ij, yA(D),ij), we can create a straight line corresponding to the ray path:(3)LA,ijs=xA,ij−s∂wA/∂x,yA,ij−s∂wA/∂y,s∂wA/∂x2+∂wA/∂y2+1x=xA,ij,y=yA,ij,LD,ijs=xD,ij−s∂wD/∂x,yD,ij−s∂wD/∂y,s∂wD/∂x2+∂wD/∂y2+1x=xD,ij,y=yD,ij,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 LB′,ij (s) and LC′,ij (s) after refraction on B′ and C′, ray positions (xB′,ij, yB′,ij) and (xC′,ij, yC′,ij) on these planes have to be established. It can be done by the calculation of LA,ij and LD,ij (x, y)-components, under the condition that z-component is equal to zB′ and zC′. This leads to the following formulas:(4)xB′,ij=LA,ijxLA,ijz=zB′=xA,ij∂wA/∂x2+∂wA/∂y2+1−zB′⋅∂wA/∂xx=xA,ij,y=xA,ij,yB′,ij=LA,ijyLA,ijz=zB′=xA,ij∂wA/∂x2+∂wA/∂y2+1−zB′⋅∂wA/∂yx=xA,ij,y=xA,ij,yC,ij=LD,ijxLD,ijz=zC′=xA,ij∂wA/∂x2+∂wA/∂y2+1−zC′⋅∂wA/∂yx=xA,ij,y=xA,ij,yC,ij=LD,ijyLD,ijz=zC′=xA,ij∂wA/∂x2+∂wA/∂y2+1−zC′⋅∂wA/∂yx=xA,ij,y=xA,ij,At each point (xB′,ij, yB′,ij) or (xC′,ij, yC′,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:(5)rB′,ij=1nrA,ij−1nrA,ijz+1−1n21−rA,ijz2⋅z^,rB′,ij=1nrD,ij−1nrD,ijz+1−1n21−rD,ijz2⋅z^,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 LB′(C′),ij will be given by the following formula:(6)LB′,ijs=xB′,ij−srB′,ijx,yB′,ij−srB′,ijy,zB′−srB′,ijz,LC′,ijs=xC′,ij−srC′,ijx,yC′,ij−srC′,ijy,zC′−srC′,ijz.The distributions of ray spots (xB,ij, yB,ij) and (xC,ij, yC,ij) can be obtained from the condition of the above ray lines intersection with B and C planes, respectively; namely, LB′C′,ijz=zBC, which gives(7)xB,ij=LB′,ijxLB′,ijz=zB=xB′,ij−rB′,ijxrB',ijzzB′−zB,xB,ij=LB′,ijyLB′,ijz=zB=yB′,ij−rB′,ijyrB′,ijzzB′−zB,xC,ij=LC′,ijxLC′,ijz=zC=xC′,ij−rC′,ijxrC′,ijzzC′−zC,yC,ij=LC′,ijyLC′,ijz=zC=yC′,ij−rC′,ijyrC′,ijzzC′−zC.These distributions define the corresponding ray spot density distributions {NB; NC}, 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:(8)IBxB,yB=NBxB,yB⋅ΔPΔSB,ICxC,yC=NCxC,yC⋅ΔPΔSC.


### 2.2. Determination of Optimal Ray-Mapping betweenB and C Planes

Irradiance transformation on the way between two planes can be equivalently understood in ray-optics sense as ray distribution change orray 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 {xC = XC(xB, yB), yC = YC(xB, yB)}, where (xB, yB) and (xC, yC) are ray coordinates at planes B and C. In the previous section, (xB,ij, yB,ij) and (xC,ij, yC,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:(9)IBxB,yBΔxBΔyB=ICXC,YCΔxCΔyC,which can also be stated as(10)IBxB,yB=ICXC,YCJXC,YC=ICXC,YC∂XC∂xB∂XC∂yB∂YC∂xB∂YC∂yB,where 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 C1-class) surface. To meet this requirement, another equation, which is called the integrability condition, has to be met: N⋅rotΝ=0, 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,(11)XC=∂Ψ∂xB;YC=∂Ψ∂yB.Considering (11) in (10), one obtains the fully nonlinear elliptic partial differential equation of Monge–Ampere type (MA PDE):(12)IBxB,yB=IC∂Ψ∂xB,∂Ψ∂yB⋅∂2Ψ∂xB2∂2Ψ∂yB2−∂2Ψ∂xB∂yB2.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 (xB,ij, yB,ij)(1) and (xC,ij, yC,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 (xB,ij, yB,ij)(1) and (xC,ij, yC,ij)(1) as follows:(13)ε1xB,yB=ICxCxB,yB,yCxB,yB⋅∂2Ψ1∂xB2∂2Ψ1∂yB2−∂2Ψ1∂xB∂yB2−IBxB,yB.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 IB1, producing slightly changed irradiance distribution at plane B, denoted by IB2. 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 IBK⟶IC transformation is “nearly” integrable. Then, by this mapping integration, the optical surface zB 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 (xB,ij, yB,ij) = (xB,ij, yB,ij)(K), i = 1, …, M, j = 1, …, N, on B plane and the corresponding collection (xC,ij, yC,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: rB,ij = [xC,ij − xB,ij, yC,ij − yB,ij, zC − zB]. Knowing the direction rB′ 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 rB′ to rB. According to the Snell law in vector notation,(14)Nij=n3⋅rB,ij−n2⋅rB′,ij,where Nij is a vector which has to be locally normal to the first freeform surface, zL1(xB, yB) at point (xB,ij, yB,ij). Having the collection of points and normals {(xB,ij, yB,ij), Nij}, one can calculate the surface shape representation: points zL1,ij and corresponding 3D point cloud (xB,ij, yB,ij, zL1,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 (xB,11, yB,11, zB) on B plane coincides with the first point of freeform surface (xB,11, yB,11, zL1,11). In other words, we state that zL1,11 = zB. It is done without losing generality, since a similar assumption could be made for any zL1,ij. The neighboring point in increasing i-index direction (xB,21, yB,21, zL1,21) can be obtained as follows. First, we define an equation of the line (ray) crossing the point (xB,21, yB,21, zB) and having the direction rB′,21:(15)l21t=xB,21+t⋅rB′,21x,yB,21+t⋅rB′,21y,zB+t⋅rB′,21z,where t is a scalar parameter corresponding to the distance from (xB,21, yB,21, zB). We can also create analytical representation of the surface (x, y, z) perpendicular to N11 at point (xB,11, yB,11, zL,11):(16)x−xB,11⋅N11x+y−yB,11⋅N11y+z−zL,11⋅N11z=0.At this point, it is possible to calculate specifict = t21, for which line l21(t21) crosses the surface defined by (16). Combining (15) and (16) leads to the following result:(17)t21=N11xxB,11−xB,21+N11yyB,11−yB,21N11⋅rB′,21.Identical formula applies for any (i, j)-index pair and both i-index and j-index increase directions:(18)ti+1j=NijxxB,ij−xB,i+1j+NijyyB,ij−yB,i+1jNij⋅rB′,i+1j,tij+1=NijxxB,ij−xB,ij+1+NijyyB,ij−yB,ij+1Nij⋅rB′,ij+1.Havingt determined, it is easy to obtain the zL1 in neighboring grid points according to the relation zL1=zB+t⋅rB′z, which specifically gives(19)zL1,i+1j=zB+NijxxB,ij−xB,i+1j+NijyyB,ij−yB,i+1jNij⋅rB′,i+1jrB′,i+1jz,zL1,ij+1=zB+NijxxB,ij−xB,ij+1+NijyyB,ij−yB,ij+1Nij⋅rB′,ij+1rB′,ij+1z.The above equations can be used to obtain iteratively the whole representation of the first freeform surfacezL1,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 sagzL2 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 surfacePC,11 has to be fixed at C plane (which translates into zL2,11 = zC). All other intersection points (PA,11, PB′,11, PB,11, PC′,11, and PD,11) are known from the previous calculations for all rays. PC,11 coordinates determine the absolute optical path length corresponding to this initial ray propagation through the system:(20)OPL11=n1PB′,11−PA,11+n2PB,11−PB′,11+n3PC,11−PB,11+n2PC′,11−PC,11+n1PD,11−PC′,11,where n1, n2, and n3 denote refractive indices of surrounding medium, glass of lenses, and medium inside beam shaper, respectively. To provide wavefront desired transformation, optical path length OPLij of any other ray (associated with ij grid point), traveling between A and D planes, should be the same as and equal to OPL11 provided in (20); additionally, the correction term resulting from input-output wavefront phase differences has to be considered, leading to(21)OPLij+wDPD,ij−wAPA,ij⋅λ2π=OPL11+wDPD,11−wAPA,11⋅λ2π.It can be further developed into the following equation:(22)n1PB′,ij−PA,ij+n2PB,ij−PB′,ij+n3PC,ij−PB,ij+n2PC′,ij−PC,ij+n1PD,ij−PC′,ij=OPL11+wDPD,11−wDPD,ij+wAPA,ij−wAPA,11⋅λ2π.In the above formula,z-coordinate of PC,ij = (xC,ij, yC,ij, zL2,ij) is the only unknown, so it can be numerically determined at desired precision. PC,ij collection creates the point cloud representing the second freeform surface zL2(xC, yC). 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.


## 2.1. Determination of Optical Field atB and C Planes

Assume that known irradianceIA(xA, yA) and wavefront wA(xA, yA) are associated with A plane. Also, the desired output at the D plane is predefined: ID(xD, yD) and wD(xD, yD). 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 IB(xB, yB) and IC(xC, yC).Irradiance is a scalar function, which describes optical power spatial distribution on a considered surface. In ray model, irradianceI(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:(1)NAxA,yA=IAxA,yA⋅ΔSAΔP,NDxD,yD=IDxD,yD⋅ΔSDΔP.The ray distribution inside elements ΔSA and ΔSD can be considered as random. As a result, we obtain the discrete distribution of points-rays (xA,ij, yA,ij) on A and (xD,ij, yD,ij) on D surface. Knowing wavefront wA(xA, yA) and wB(xB, yB), local ray directions on both planes can be represented as unit normal vectors to those wavefronts:(2)rA,ij=−∂wA/∂x,−∂wA/∂y,1∂wA/∂x2+∂wA/∂y2+1x=xA,ij,y=yA,ij,rD,ij=−∂wD/∂x,−∂wD/∂y,1∂wD/∂x2+∂wD/∂y2+1x=xD,ij,y=yD,ij.Now, at any ray spot (xA(D),ij, yA(D),ij), we can create a straight line corresponding to the ray path:(3)LA,ijs=xA,ij−s∂wA/∂x,yA,ij−s∂wA/∂y,s∂wA/∂x2+∂wA/∂y2+1x=xA,ij,y=yA,ij,LD,ijs=xD,ij−s∂wD/∂x,yD,ij−s∂wD/∂y,s∂wD/∂x2+∂wD/∂y2+1x=xD,ij,y=yD,ij,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 LB′,ij (s) and LC′,ij (s) after refraction on B′ and C′, ray positions (xB′,ij, yB′,ij) and (xC′,ij, yC′,ij) on these planes have to be established. It can be done by the calculation of LA,ij and LD,ij (x, y)-components, under the condition that z-component is equal to zB′ and zC′. This leads to the following formulas:(4)xB′,ij=LA,ijxLA,ijz=zB′=xA,ij∂wA/∂x2+∂wA/∂y2+1−zB′⋅∂wA/∂xx=xA,ij,y=xA,ij,yB′,ij=LA,ijyLA,ijz=zB′=xA,ij∂wA/∂x2+∂wA/∂y2+1−zB′⋅∂wA/∂yx=xA,ij,y=xA,ij,yC,ij=LD,ijxLD,ijz=zC′=xA,ij∂wA/∂x2+∂wA/∂y2+1−zC′⋅∂wA/∂yx=xA,ij,y=xA,ij,yC,ij=LD,ijyLD,ijz=zC′=xA,ij∂wA/∂x2+∂wA/∂y2+1−zC′⋅∂wA/∂yx=xA,ij,y=xA,ij,At each point (xB′,ij, yB′,ij) or (xC′,ij, yC′,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:(5)rB′,ij=1nrA,ij−1nrA,ijz+1−1n21−rA,ijz2⋅z^,rB′,ij=1nrD,ij−1nrD,ijz+1−1n21−rD,ijz2⋅z^,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 LB′(C′),ij will be given by the following formula:(6)LB′,ijs=xB′,ij−srB′,ijx,yB′,ij−srB′,ijy,zB′−srB′,ijz,LC′,ijs=xC′,ij−srC′,ijx,yC′,ij−srC′,ijy,zC′−srC′,ijz.The distributions of ray spots (xB,ij, yB,ij) and (xC,ij, yC,ij) can be obtained from the condition of the above ray lines intersection with B and C planes, respectively; namely, LB′C′,ijz=zBC, which gives(7)xB,ij=LB′,ijxLB′,ijz=zB=xB′,ij−rB′,ijxrB',ijzzB′−zB,xB,ij=LB′,ijyLB′,ijz=zB=yB′,ij−rB′,ijyrB′,ijzzB′−zB,xC,ij=LC′,ijxLC′,ijz=zC=xC′,ij−rC′,ijxrC′,ijzzC′−zC,yC,ij=LC′,ijyLC′,ijz=zC=yC′,ij−rC′,ijyrC′,ijzzC′−zC.These distributions define the corresponding ray spot density distributions {NB; NC}, 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:(8)IBxB,yB=NBxB,yB⋅ΔPΔSB,ICxC,yC=NCxC,yC⋅ΔPΔSC.


## 2.2. Determination of Optimal Ray-Mapping betweenB and C Planes

Irradiance transformation on the way between two planes can be equivalently understood in ray-optics sense as ray distribution change orray 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 {xC = XC(xB, yB), yC = YC(xB, yB)}, where (xB, yB) and (xC, yC) are ray coordinates at planes B and C. In the previous section, (xB,ij, yB,ij) and (xC,ij, yC,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:(9)IBxB,yBΔxBΔyB=ICXC,YCΔxCΔyC,which can also be stated as(10)IBxB,yB=ICXC,YCJXC,YC=ICXC,YC∂XC∂xB∂XC∂yB∂YC∂xB∂YC∂yB,where 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 C1-class) surface. To meet this requirement, another equation, which is called the integrability condition, has to be met: N⋅rotΝ=0, 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,(11)XC=∂Ψ∂xB;YC=∂Ψ∂yB.Considering (11) in (10), one obtains the fully nonlinear elliptic partial differential equation of Monge–Ampere type (MA PDE):(12)IBxB,yB=IC∂Ψ∂xB,∂Ψ∂yB⋅∂2Ψ∂xB2∂2Ψ∂yB2−∂2Ψ∂xB∂yB2.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 (xB,ij, yB,ij)(1) and (xC,ij, yC,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 (xB,ij, yB,ij)(1) and (xC,ij, yC,ij)(1) as follows:(13)ε1xB,yB=ICxCxB,yB,yCxB,yB⋅∂2Ψ1∂xB2∂2Ψ1∂yB2−∂2Ψ1∂xB∂yB2−IBxB,yB.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 IB1, producing slightly changed irradiance distribution at plane B, denoted by IB2. 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 IBK⟶IC transformation is “nearly” integrable. Then, by this mapping integration, the optical surface zB 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 (xB,ij, yB,ij) = (xB,ij, yB,ij)(K), i = 1, …, M, j = 1, …, N, on B plane and the corresponding collection (xC,ij, yC,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: rB,ij = [xC,ij − xB,ij, yC,ij − yB,ij, zC − zB]. Knowing the direction rB′ 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 rB′ to rB. According to the Snell law in vector notation,(14)Nij=n3⋅rB,ij−n2⋅rB′,ij,where Nij is a vector which has to be locally normal to the first freeform surface, zL1(xB, yB) at point (xB,ij, yB,ij). Having the collection of points and normals {(xB,ij, yB,ij), Nij}, one can calculate the surface shape representation: points zL1,ij and corresponding 3D point cloud (xB,ij, yB,ij, zL1,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 (xB,11, yB,11, zB) on B plane coincides with the first point of freeform surface (xB,11, yB,11, zL1,11). In other words, we state that zL1,11 = zB. It is done without losing generality, since a similar assumption could be made for any zL1,ij. The neighboring point in increasing i-index direction (xB,21, yB,21, zL1,21) can be obtained as follows. First, we define an equation of the line (ray) crossing the point (xB,21, yB,21, zB) and having the direction rB′,21:(15)l21t=xB,21+t⋅rB′,21x,yB,21+t⋅rB′,21y,zB+t⋅rB′,21z,where t is a scalar parameter corresponding to the distance from (xB,21, yB,21, zB). We can also create analytical representation of the surface (x, y, z) perpendicular to N11 at point (xB,11, yB,11, zL,11):(16)x−xB,11⋅N11x+y−yB,11⋅N11y+z−zL,11⋅N11z=0.At this point, it is possible to calculate specifict = t21, for which line l21(t21) crosses the surface defined by (16). Combining (15) and (16) leads to the following result:(17)t21=N11xxB,11−xB,21+N11yyB,11−yB,21N11⋅rB′,21.Identical formula applies for any (i, j)-index pair and both i-index and j-index increase directions:(18)ti+1j=NijxxB,ij−xB,i+1j+NijyyB,ij−yB,i+1jNij⋅rB′,i+1j,tij+1=NijxxB,ij−xB,ij+1+NijyyB,ij−yB,ij+1Nij⋅rB′,ij+1.Havingt determined, it is easy to obtain the zL1 in neighboring grid points according to the relation zL1=zB+t⋅rB′z, which specifically gives(19)zL1,i+1j=zB+NijxxB,ij−xB,i+1j+NijyyB,ij−yB,i+1jNij⋅rB′,i+1jrB′,i+1jz,zL1,ij+1=zB+NijxxB,ij−xB,ij+1+NijyyB,ij−yB,ij+1Nij⋅rB′,ij+1rB′,ij+1z.The above equations can be used to obtain iteratively the whole representation of the first freeform surfacezL1,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 sagzL2 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 surfacePC,11 has to be fixed at C plane (which translates into zL2,11 = zC). All other intersection points (PA,11, PB′,11, PB,11, PC′,11, and PD,11) are known from the previous calculations for all rays. PC,11 coordinates determine the absolute optical path length corresponding to this initial ray propagation through the system:(20)OPL11=n1PB′,11−PA,11+n2PB,11−PB′,11+n3PC,11−PB,11+n2PC′,11−PC,11+n1PD,11−PC′,11,where n1, n2, and n3 denote refractive indices of surrounding medium, glass of lenses, and medium inside beam shaper, respectively. To provide wavefront desired transformation, optical path length OPLij of any other ray (associated with ij grid point), traveling between A and D planes, should be the same as and equal to OPL11 provided in (20); additionally, the correction term resulting from input-output wavefront phase differences has to be considered, leading to(21)OPLij+wDPD,ij−wAPA,ij⋅λ2π=OPL11+wDPD,11−wAPA,11⋅λ2π.It can be further developed into the following equation:(22)n1PB′,ij−PA,ij+n2PB,ij−PB′,ij+n3PC,ij−PB,ij+n2PC′,ij−PC,ij+n1PD,ij−PC′,ij=OPL11+wDPD,11−wDPD,ij+wAPA,ij−wAPA,11⋅λ2π.In the above formula,z-coordinate of PC,ij = (xC,ij, yC,ij, zL2,ij) is the only unknown, so it can be numerically determined at desired precision. PC,ij collection creates the point cloud representing the second freeform surface zL2(xC, yC). 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:zA = 0 mm, zB′ = 25 mm, zB = 30 mm, zC = 135 mm, zC′ = 160 mm, zD = 500 mm, n1 = n3 = 1 (air is both the surrounding medium and the medium between lenses), n2 = 1.51. The described methodology was implemented in Matlab to calculate the numerical representations of both freeform surfaces (Figure 3).Figure 3
Contour plots (up) and CAD renders (down) showing the calculated freeform surfaces sag data:zL1 (a) and zL2 (b).
(a)(b)(c)(d)Additionally, wavefront error at planeD 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:(23)OPDij=PB′,ij−PA,ij+n2PB,ij−PB′,ij+PC,ij−PB,ij+n2PC′,ij−PC,ij+PD,ij−PC′,ij.It appeared that wavefront flatness was obtained at very high precision; OPD did not exceed 10−12 mm (Figure 4).Figure 4
Optical path difference (wavefront error) of the designed beam shaper.In order to verify the designed beam shaper performance, the data was transferred asgrid 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.Figure 5
Screenshot from OpticStudio experiment.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.Figure 6
Spot diagram representing ray distribution inD plane (z = 500 mm).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 (Figure7).Figure 7
Spot diagrams representing ray distributions beyondD plane: z = 10 m (a), z = 100 m (b).
(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 distancezr: xrd, yrd and its real position xr, yr, the ray direction error δ (vector) can be calculated from the following formula:(24)δ=1zxrd−xr,yrd−yr.Since a ray is always perpendicular to a wavefront Φ, ray direction error determines wavefront slope (wavefront gradient):(25)∇ϕ=δ.Taking into account (24) and (25), it can be stated that(26)ϕx=xr−xrdz,ϕy=yr−yrdz,where Φ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.Figure 8
Fried geometry for thei, j ray. Wavefront gradient components Φx,ij and Φy,ij are obtained from (26) and serve to calculate wavefront phase Φij distribution.Assuming perfect ray positions, which form a triangular shape defining a border of beam footprint, wavefront at half grid points inside can be approximated as(27)ϕx,ij=12ϕi+1/2,j−1/2−ϕi−1/2,j−1/2+ϕi+1/2,j+1/2−ϕi−1/2,j+1/2dx,(28)ϕy,ij=12ϕi−1/2,j+1/2−ϕi−1/2,j−1/2+ϕi+1/2,j+1/2−ϕi+1/2,j−1/2dy,where dx and dy 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 N2 ×1 Φ-vector. In the same way, N×N×2 array of known wavefront slopes [Φx; Φy] can be stacked into N2×2 g-array. Then (27) and (28) can be transformed into matrix form:(29)g=Γϕ,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:(30)Γϕ−g22.Following [29], the minimum-norm solution is given by(31)ϕ^=ΓTΓ†ΓTg,where superscript T corresponds to matrix transposition and † denotes its pseudoinverse. To apply the discussed method, ray positions at z = 100 m were considered as xr, yr and desired ray positions at D plane as xrd,yrd. 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).Figure 9
Wavefront error (lenses resolved as 600 × 600 grids).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 (Figure10); 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).Figure 10
Wavefront error for increased resolution (lenses resolved as 1000 × 1000 grids).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.



---

*Source: 2896593-2020-08-01.xml*