## Abstract

Many laser applications require specific irradiance distributions to ensure optimal performance. Geometric optical design methods based on numerical calculation of two plano-aspheric lenses have been thoroughly studied in the past. In this work, we present an alternative new design approach based on functional differential equations that allows direct calculation of the rotational symmetric lens profiles described by two-point Taylor polynomials. The formalism is used to design a Gaussian to flat-top irradiance beam shaping system but also to generate a more complex dark-hollow Gaussian (donut-like) irradiance distribution with zero intensity in the on-axis region. The presented ray tracing results confirm the high accuracy of both calculated solutions and emphasize the potential of this design approach for refractive beam shaping applications.

© 2014 Optical Society of America

## 1. Introduction

Many commercial, medical and scientific applications of the laser have been developed since its invention. Materials processing applications, like cutting or welding; and medical applications, like corneal surgery are just a few examples. Many of these applications require a specific beam irradiance distribution to ensure an optimal performance [1]. There are several optical technologies that can be used for laser beam shaping, including simple apertures, diffractive optical elements [2], lenses [3] and mirrors [4]. The right choice depends on the optical design problem. Refractive laser beam shapers are commonly used due to their high optical efficiency and simple structure [5].

Often, it is possible to apply geometrical methods to shape a laser beam profile. This common design approach is based on the ray mapping between the input and the output plane which takes into account conservation of energy, a constant optical path length of the geometric wave-front between the input and output planes, and the ray trace equations. Geometric ray mapping-based designs with two plano-aspheric lenses have been thoroughly studied and improved [6–10]. Figure 1 shows a typical Galilean configuration (without internal focusing) of a refractive laser beam shaper consisting of two plano-aspheric lenses transforming a Gaussian input into a flattop output intensity distribution with a gradual roll-off from the uniform to null region. The indicated geometric ray mapping then describes where an initial parallel ray with radial distance *r* is refracted at the first aspherical surface to be redirected in parallel direction with radial distance *R* at the second aspherical surface.

Even though analytic expressions of ray mapping functions do exist for some specific cases, the surface profiles of the lenses are still calculated numerically. The main idea of most common geometrical optics design approaches is to describe the laser beam shaping design problem by mathematical equations and then use numerical techniques to solve the particular equation. Solutions have been found in the form of an integral equation, by Kreuzer, ordinary differential equation, by Rhodes and Shealy, and Monge-Ampère equation, by Oliker. Usually, a set of discrete data points is obtained using numerical calculation techniques and the final profiles of the lenses are constructed by either polynomial or spline fitting.

In this work, we present a novel different design approach that describes the entire optical design problem based on functional differential equations. This design approach is related to an analytic design method for optimal imaging that enables the coupling of three ray-bundles with only two free-form lens surfaces in two and three dimensions [11, 12]. However, the original design approach worked so far only for spherical and plane wave-fronts. Therefore, further extensions are necessary and key to properly incorporate the ray mapping function.

Starting from the example of a Gaussian to flat-top irradiance beam shaping system in Fig. 1, a set of two functional differential equations are derived in Sec. 2. The mathematical problem is then converted into a set of linear algebraic equations that allows calculating analytic symbolic expressions for the lens profiles *f* and *g* described by two-point Taylor series expansions. To further demonstrate the versatility of this new approach, a second design problem is solved in Sec. 3 to generate a more complex dark-hollow Gaussian (donut-like) irradiance profile with zero intensity in the on-axis region. Here, the mapping function is given by a fractional Taylor series which has direct impact on the surface representation as well. The presented ray tracing evaluation confirm the high accuracy of both calculated solutions and emphasize the potential of this design approach for refractive beam shaping applications.

## 2. Analytic solution to convert a Gaussian into a flat-top irradiance distribution

If the lenses are rotationally symmetric and described by surfaces of revolution, this means that the problem can be examined in a plane. The input irradiance distribution is given by the Gaussian function

*e*

^{2}width

*w*

_{0}. To avoid diffraction effects in regions of abrupt irradiance change and to increase the distance over which the uniform irradiance distribution can be used [13], the output intensity is chosen as a flattened Lorentzian (FL) distribution which has a flat central region and a gradual roll-off to the null region, shown in Fig. 1. Where

*q*denotes the beam shape parameter and

*R*is the beam width parameter of the flattened Lorentzian [13]. The input encircled energy is calculated as a function of

_{FL}*r*as

*R*as

*R*as a function of the given input coordinate

*r*such that the encircled energies

*A*and

*B*are equal, that is

*ε*is equal to +1 for a Galilean and −1 for a Keplerian configuration. An explicit mapping function

*R*(

*r*) that transforms a Gaussian into a flattened Lorentzian irradiance distribution is shown as dotted line in Fig. 2(a) for

*ε*= 1,

*w*

_{0}= 2.366mm,

*R*= 3.25mm and

_{FL}*q*= 16.

Now that the output coordinate value *R* can be calculated for every input coordinate *r*, it is possible to derive an integral representation for the sag of each aspherical lens that can be solved numerically following the work of Kreuzer, Shealy and Hoffnagle [13].

Instead of following this well-known design approach, an alternative solution strategy will be explained hereafter. As both incoming and outgoing geometric wave-fronts are plane (i.e. all rays are parallel), the two plane surfaces in Fig. 1 have no optical power and can be neglected by assuming that the in- and output rays are immersed in a medium with refractive index *n*_{2}. This is shown in Fig. 2(b). Consider an arbitrary light ray path. The ray is emitted from a point *v⃗*_{0} on the plane incident wave-front and intersects the first surface at a point *p⃗* = (*r*, *f*(*r*)), where it is refracted towards (*R*(*r*), *g*(*R*(*r*))) on the second surface and refracted again towards the plane outgoing wave-front (parallel to the *z*-axis). The path of the light rays are governed by Fermat’s principle which states that the optical path length (OPL) between two points is an extremum along the light ray containing these points. The OPL between points *v⃗*_{0}, *p⃗* and (*R*(*r*), *g*(*R*(*r*))) is given by

*n⃗*

_{0}denotes the directional vector of the wave-front in

*z*-direction. If the points

*v⃗*

_{0}and (

*R*(

*r*),

*g*(

*R*(

*r*))) are fixed, the only free parameter that can vary to achieve an extremum for

*V*

_{1}is the point

*p⃗*and hence

*r*on the first lens surface. Fermat’s principle thus implies that where the partial derivative indicates that

*R*is held fixed. Furthermore, the slope of the first aspheric surface at point (

*r*,

*f*(

*r*)) must be equal to the slope on the second aspheric surface at point (

*R*,

*g*(

*R*)) as a result of symmetry of the input and output wave-fronts. This can be expressed in a second condition where the partial derivative

*∂*(

_{r}g*r*) is evaluated at position

*R*(

*r*). This specific formulation highlights the functional character of the optical design problem. As the mapping function

*R*is a function of

*r*and appears as an argument in the unknown function

*g*, Eqs. (7) and (8) form a pair of functional differential equations for two unknown functions

*f*(

*r*) and

*g*(

*R*(

*r*)).

Suppose (*f*, *g*) is an analytic and smooth solution to the functional differential equations (7) and (8), Taylor’s theorem implies that the functions must be infinitely differentiable and have a power-series representation. In our previous work [11], a Taylor series approach has been successfully used to solve a related set of functional differential equations by solving equations

*n*

^{th}order, evaluated at that point

*r*

_{1}. If the mapping function is now reconstructed from its (exact) derivatives at

*r*

_{1}, the Taylor series of this mapping function would not converge for all values

*r*(see Fig. 2(a)) and the same applies for the lens profiles

*f*and

*g*. Therefore, as alternative, a Taylor series expansion about two points

*r*

_{1}and

*r*

_{2}can be defined as

*k*+ 1)

^{th}degree polynomial approximation of the mapping function

*R*(

*r*), if the coefficients

*a*and

_{i}*b*do fulfill that

_{i}*∂*is equal to

^{i}T/∂r^{i}*∂*at

^{i}R/∂r^{i}*r*

_{1}and

*r*

_{2}for

*i*= 0, 1,..,

*k*[14]. The two-point Taylor series approximation of the mapping function converges well for all values

*r*and is added to the graph in Fig. 2(a). The selection of points

*r*

_{1}and

*r*

_{2}depends strongly on the shape of the considered mapping function. To find a suitable pair of points, it has been found sufficient to verify that the two-point Taylor series approximation of the mapping function converges well for all values

*r*. This was always the case even for manually selected

*r*

_{1}values close to the optical axis and

*r*

_{2}values close to the maximum radius. There is no strict mathematical proof that the same convergence also applies for the lens profiles. However, all evaluated designs demonstrated this behavior so far. The results presented in this work will provide strong evidence that the solutions for the lens profiles are analytic and smooth as well. Therefore, the surface profile functions can be expressed by two-point Taylor series expansions

*r*

_{1},

*z*

_{1}) and (

*r*

_{2},

*z*

_{2}), and (

*r*

_{3},

*z*

_{3}) and (

*r*

_{4},

*z*

_{4}) respectively. The two values

*r*

_{3}=

*R*(

*r*

_{1}) and

*r*

_{4}=

*R*(

*r*

_{2}) are given by the mapping function

*R*(

*r*), respectively. The ray paths of both initial rays are shown in Fig. 2(b). For the first ray,

*z*

_{1}can be freely chosen without loss of generality, whereas the value

*z*

_{3}fixes the distance between the two lens surfaces. For the second ray, the position

*z*

_{2}can not be determined analytically and represents a single parameter that links the two local solutions at

*r*

_{1}and

*r*

_{2}. This fact will be explained in more detail at the beginning of Sec. 2.1. The position

*z*

_{4}then follows from a constant OPL. All four slope values

*m*

_{1},

*m*

_{2},

*m*

_{3}and

*m*

_{4}at points

*r*(

_{i}*i*= 1..4) follow from the ray mapping relationship and Snell’s law.

The initial conditions

*D*= 0 for

_{i}*i*= 1, 2 and provide general solutions for the initial two-point Taylor series coefficients depending upon the single variable

*z*

_{2}. In ascending order, it is now possible to calculate (2

*k*+ 1)

^{th}order Taylor polynomials of

*f*(

*r*) and

*g*(

*r*) by calculating the derivatives and evaluating equations

*r*

_{1}and

*r*

_{2}. For the case

*n*= 0, the initial Taylor series coefficients

*p*

_{0},

*p*

_{1},

*q*

_{0},

*q*

_{1},

*u*

_{0},

*u*

_{1},

*v*

_{0}and

*v*

_{1}as defined in Eqs. (13) solve Eqs. (14).

For every higher order *n* ≥ 1, the set of equations (14) results in four linear algebraic equations for two-point Taylor series coefficients *p _{i}*,

*q*,

_{i}*u*and

_{i}*v*. By sorting and combining all terms, the equations can be expressed as a compact matrix equation for particular two-point Taylor series coefficients for arbitrary

_{i}*n*≥ 1. The elements of the linear algebraic equations consist of mathematical expressions only dependent on previously calculated Taylor series coefficients. Therefore, all two-point Taylor series coefficients can be calculated as exact symbolic expressions in ascending order - step by step. So far, no approximations have been made. The general solution scheme for

*n*≥ 1 allows to calculate successive Taylor series coefficients of

*f*(

*r*) and

*g*(

*r*) up to a very high order. The inevitable truncation of the infinite sum of terms will be the only approximation made. The radii of convergence for the Taylor expansions

*f*(

*r*) and

*g*(

*r*) are very important, as they indicate the maximum aperture that can be achieved for a given set of initial values. For the examples considered in this work, the radii of convergence are larger than the radial coordinate 0 <

*r*<

*r*of the lenses.

_{max}#### 2.1. Evaluation of the analytic solution

The parameters of the in- and output intensity distributions are *w*_{0} = 2.366mm, *R _{FL}* = 3.25mm and

*q*= 16 as before. The two points of the Taylor series are

*r*

_{1}= 0.2mm and

*r*

_{2}= 3.6mm. The chosen distance between

*z*

_{1}= 0mm and

*z*

_{3}= 16mm defines the lens separation. For a chosen refractive index

*n*

_{2}= 1.5, this initial ray path construction determines the OPL. After applying the constant OPL condition for point

*r*

_{2}, there still exists an infinite number of solutions that solve the design problem locally in the neighborhood of

*r*

_{2}. In order to link the two local Taylor expansion solutions at

*r*

_{1}and

*r*

_{2}, the correct relative positioning in

*z*-direction of the two solution parts has to be identified. It is very important to note that both local solutions at the respective points are exact and calculated symbolically, only the relative position of the two parts of the overall solution has to be adjusted numerically. To find the appropriate value of

*z*

_{2}, the implemented calculation of the Taylor coefficients for the lens profiles is defined as a function in

*Mathematica*. After calculation of the Taylor series coefficients up to a certain order

*k*, 10 equidistant rays between 0 mm and 4 mm (corresponds to ≈ 99.7% encircled input energy) are traced and evaluated. The spatial deviation of each ray is then defined by the mapping error

*R′*−

_{i}*R*(

*r*). Where

_{i}*R′*

*is the calculated intersection of the*

_{i}*i*

^{th}ray with the surface profile

*g*and

*R*(

*r*) is the ideal position defined by the mapping function. As a second figure of merit, the variation in optical path length is also calculated for all rays. As one possible overall merit function,

_{i}*m*(

*z*

_{2}) is given by the sum of the two root mean square values of the spatial deviation and variation in OPL of the 10 rays over the full entrance aperture.

Figure 3 shows the evaluation of the merit function *m*(*z*_{2}) as a function of variable *z*_{2} for values between −0.5 and 0 mm and for a step size of 0.0025 mm. A logarithmic plot is chosen to better cover the wide range of values.

The curve in parameter space is characterized by a convex curvature which makes it easy and fast to find the optimal relative position of *z*_{2} using a local optimization in *Mathematica*.

To quantitatively evaluate the accuracy of the derived solution, the spatial and angular deviation is calculated for 200 rays between 0 and 4 mm along *r*-direction. This evaluation is performed for 7^{th} (k=3), 13^{th} (k=6) and 19^{th} (k=9) order Taylor polynomials for *f* and *g*, respectively. Figure 4(a) shows the variation in optical path length (OPL) and Fig. 4(b) shows the spatial mapping error for 3 calculated solutions of different order. The root mean square values of the variation in OPL are 1.647μm (k=3), 0.617μm (k=6) and 0.058μm (k=9). The root mean square values of the mapping error are 0.559 × 10^{−2} (k=3), 0.061 × 10^{−2} (k=6) and 0.009 × 10^{−2} (k=9), respectively. For both graphs, the variation in OPL and mapping error reduce very well with an increasing order underlining the convergence and achievable high accuracy of the derived solution.

Furthermore, the maximum collimation error after refraction at the second plane surface is evaluated for the 3 different solutions. The maximum collimation error is 2.167mrad (k=3), 1.124mrad (k=6) and 0.265mrad (k=9). The achieved collimation quality for k=9 is clearly better than the maximum collimation error of 0.8mrad in [7]. Further improvements, if necessary, can be achieved by increasing the order of the Taylor polynomials and/or by adjusting the positions of points *r*_{1} and *r*_{2}, respectively.

## 3. Analytic solution for Gaussian to dark-hollow Gaussian refractive beam shapers

The existence of a closed-form expression for the mapping function was important for the design approach in Sec. 2 to work. The following example will demonstrate that this is not a necessary condition. The optical design problem of a refractive beam shaper that transforms a Gaussian into a dark-hollow Gaussian (DHG) intensity profile with zero intensity at the optical axis does not allow to calculate a closed-form expression for the mapping function, as will be shown in the following. Focused small ring patterns are of interest for some industrial applications such as laser cladding technology. Figure 5 shows the cross-section of the correspondent setup together with the in- and output intensity profiles.

The input irradiance distribution is again given by the Gaussian function as defined in Eq. (1). The output intensity distribution is given by

*n*denotes the order of the DHG beam and

*H*

_{0}is a constant that is adjusted to conserve the overall energy. For

*n*= 0, Eq. (15) reduces to the fundamental beam with beam waist

*w*

_{1}[15].

As done in the previous Section, it is possible to determine the relationship between output radial distance *R* and input coordinate *r* such that the encircled energy *B*(*R*) at the output is equal to the encircled energy *A*(*r*) at the input. The input encircled energy was given in Eq. (3) and the output encircled energy is

*R*(

*r*) can not be solved as a closed-form expression. However, it is possible to derive a closed-form expression for the inverse mapping function

*r*(

*R*) =

*R*

^{−1}(

*r*). The mapping function

*R*(

*r*) is then given by applying Lagrange inversion theorem that provides the Taylor series expansion of the inverse function of an analytic function. Both mapping function

*r*(

*R*) and its inverse function

*R*(

*r*) are shown in Fig. 6(a); for parameters

*w*

_{0}= 2.366mm,

*w*

_{1}= 2mm and

*n*= 1.

Interestingly, the inverse mapping function *R*(*r*) is then described by a Taylor series expansion

*R′*(

*r*) for

*r*→ 0 diverges, as the function intersects the horizontal axis perpendicular at the origin (see dotted line in Fig. 6(a)). When trying to use the same solution scheme for the approximated mapping function

*R*(

*r*) as in the previous section, the problem appears that even for a two-point Taylor series representation of the surfaces the profiles will still not converge for all values of

*r*. This problem is directly linked to the fractional exponents of the mapping function

*R*(

*r*), which will become more clear later.

To solve this design problem, it is better to solve its inverse problem: where a dark-hollow Gaussian is transformed into a Gaussian beam. This inverse optical design problem is identical due to the symmetry of the optical design where both in- and outgoing geometric wave-fronts are plane. The inverted design is shown in Fig. 6(b), where the reversed directions of the rays are indicated. As done previously, this time from right to left, the OPL between points *v⃗*_{0}, *p⃗* = (*R*, *g*(*R*)) and (*r*(*R*), *f*(*r*(*R*))) is given by

*r*is held fixed. To ensure recollimation of the rays, the second condition is given by where the partial derivative

*∂*

_{R}*f*(

*R*) is evaluated at position

*r*(

*R*). The surface

*g*(

*R*) is defined by a two-point Taylor series expansion of (2

*k*+ 1)

^{th}order

*R*

_{1}and

*R*

_{2}in ascending order

*g*(

*R*) in ascending order up to high orders. However, a similar two-point Taylor polynomial of function

*f*still does not converge for all points

*R*. From Eq. (21) follows that it is also possible to calculate the second surface profile

*f*(

*r*(

*R*)) from the integral equation for known functions

*g*(

*R*) and

*R*(

*r*). By substituting the mapping function

*R*(

*r*) from Eq. (18) as argument in the first derivative

*g′*(

*R*), the unknown function

*f*(

*R*) can be exactly calculated via symbolic integration. The function

*f*(

*R*) is then described by a polynomial function

*f*

_{0}is calculated by applying the constant OPL condition at point

*R*

_{2}. Interestingly, the coefficient

*f*

_{1}is almost zero (less than 10

^{−8}) and identical to the slope value

*g′*(0), providing smooth lens shapes in the central on-axis region. Due to the substitution of the mapping function

*R*(

*r*) in

*g′*(

*R*),

*f*(

*R*) easily results in a very high order polynomial. Even though it is not necessary, this function can be truncated to a moderate order in the range of function

*g*(

*R*) without changing the performance.

#### 3.1. Evaluation of the analytic solution

The performance evaluation of the analytic solution for Gaussian to dark-hollow Gaussian (DHG) beam shaping is carried out analogously to Sec. 2.1. The parameters of the in- and output intensity distributions are *w*_{0} = 2.366mm, *w*_{1} = 2mm and *n* = 1. The two manually selected points of the Taylor series are *R*_{1} = 0.3mm and *R*_{2} = 4.2mm. By fixing the distance between *Z*_{1} = 0mm and *Z*_{2} = 16mm and the refractive index *n*_{2} = 1.5, the remaining parameter *Z*_{2} that adjusts the relative positioning in *z*-direction is optimized for an identical merit function as defined in Sec. 2.1. Figure 7(a) shows the variation in OPL and Fig. 7(b) shows the spatial mapping error for 3 different solutions of different order. All calculations have been repeated for 200 rays between 0 and 4 mm along *r*-direction (which corresponds to ≈ 99.7% encircled input energy). The evaluation is performed for 7^{th} (k=3), 11^{th} (k=5) and 15^{th} (k=7) order Taylor polynomials for *g*, and *f* has been truncated to 21^{th} order.

The root mean square values of the variation in OPL are 0.053μm (k=3), 0.064μm (k=5) and 0.025μm (k=7). The root mean square values of the mapping error are 0.792 × 10^{−2} (k=3), 0.776 × 10^{−2} (k=5) and 0.164 × 10^{−2} (k=7), respectively. As in the previous design example, the variation in OPL and the mapping error reduce very well with an increasing Taylor series order, confirming the convergence and high accuracy of the derived analytic solution. Furthermore, the maximum collimation error after refraction at the second plane surface is evaluated for the 3 different solutions. The maximum collimation error is 0.744mrad (k=3), 0.063mrad (k=5) and 0.058mrad (k=7). If required, an even higher accuracy can be achieved by increasing the order of the Taylor polynomials and/or by fine-tuning the radial positions *R*_{1} and *R*_{2}.

## 4. Conclusion

Within the scope of this work, it has been shown how the optical design of refractive laser beam shaping systems based on two plano-aspheric lenses can be described by means of a new mathematical formulation using functional differential equations. The derived solution scheme has been successfully verified for two specific design tasks and allows an increasingly accurate calculation of the lens profiles described by high order Taylor polynomials. In case of the flattened Lorentzian beam shaping system, a maximum collimation error value of 0.265mrad has been achieved. If required, a higher accuracy seems feasible with this new design approach.

In case of the dark-hollow Gaussian output intensity distribution, a maximum collimation error value of 0.058mrad has been achieved. Furthermore, the design method established a direct link between the mathematical structure of the mapping function and the mathematical function describing the first optical surface (fractional exponents). This deeper insight provided by the optical design approach could prove beneficial when it comes to designing and manufacturing of such optical systems. The description of the lens surfaces with analytic functions allows imposing boundary conditions (e.g., for the local curvature of the surfaces) directly in the design process to meet manufacturing requirements.

Particularly noteworthy is the successful implementation of a two-point Taylor series expansion that could be also of interest for other optical designs, for example in imaging systems [16]. Here, a multi-point Taylor series approach could be used to generalize the optical design for more than two surfaces that allows the simultaneous coupling of more than three ray-bundles.

A related functional differential equation based optical design approach has been already applied successfully to design a system with movable free-form optics for a new class of solar concentrators [17]. This already demonstrated possibility to incorporate movement in the mathematical description could be useful for the design of dynamic refractive laser beam shaping systems. In such a system, the relative motion of the lenses will be used to dynamically change the shape or profile of the laser beam: a promising functionality for the laser materials processing industry. This dynamic optical design concept will be addressed in future work.

## Acknowledgments

The work reported in this paper is supported in part by the Research Foundation - Flanders (FWO-Vlaanderen) that provides a post-doctoral grant to Fabian Duerr, and in part by the IAP-BELSPO grant IAP P7-35 photonics@be, the Industrial Research Funding (IOF), Methusalem, and the OZR of the Vrije Universiteit Brussel.

## References and links

**1. **F. M. Dickey, S. C. Holswade, and D. L. Shealy, *Laser Beam Shaping Applications* (CRC Press, 2005). [CrossRef]

**2. **D. Palima and J. Glückstad, “Gaussian to uniform intensity shaper based on generalized phase contrast,” Opt. Express **16**, 1507–1516 (2008). [CrossRef] [PubMed]

**3. **J. A. Hoffnagle and C. M. Jefferson, “Design and performance of a refractive optical system that converts a Gaussian to a flattop beam,” Appl. Opt. **39**, 5488–5499 (2000). [CrossRef]

**4. **V. Oliker, “Optical design of freeform two-mirror beam-shaping systems,” J. Opt. Soc. Am. A **24**, 3741–3752 (2007). [CrossRef]

**5. **L. Romero and F. Dickey, “Lossless laser beam shaping,” J. Opt. Soc. Am. A **13**, 751–760 (1996). [CrossRef]

**6. **P. W. Rhodes and D. L. Shealy, “Refractive optical systems for irradiance redistribution of collimated radiation: their design and analysis,” Appl. Opt. **19**, 3545–3553 (1980). [CrossRef] [PubMed]

**7. **J. A. Hoffnagle and C. M. Jefferson, “Beam shaping with a plano-aspheric lens pair,” Opt. Eng. **42**, 3090–3099 (2003). [CrossRef]

**8. **H. Ma, Z. Liu, P. Jiang, X. Xu, and S. Du, “Improvement of Galilean refractive beam shaping system for accurately generating near-diffraction-limited flattop beam with arbitrary beam size,” Opt. Express **19**, 13105–13117 (2011). [CrossRef] [PubMed]

**9. **J. Rubinstein and G. Wolansky, “Intensity control with a free-form lens,” J. Opt. Soc. Am. A **24**, 463–469 (2007). [CrossRef]

**10. **Z. Feng, L. Huang, G. Jin, and M. Gong, “Designing double freeform optical surfaces for controlling both irradiance and wavefront,” Opt. Express **21**, 28693–28701 (2013). [CrossRef]

**11. **F. Duerr, P. Benítez, J. C. Miñano, Y. Meuret, and H. Thienpont, “Analytic design method for optimal imaging: coupling three ray sets using two free-form lens profiles,” Opt. Express **20**, 5576–5585 (2012). [CrossRef] [PubMed]

**12. **F. Duerr, P. Benítez, J. C. Miñano, Y. Meuret, and H. Thienpont, “Analytic free-form lens design in 3D: coupling three ray sets using two lens surfaces,” Opt. Express **20**, 10839–10846 (2012). [CrossRef] [PubMed]

**13. **D. L. Shealy and J. A. Hoffnagle, “Review: design and analysis of plano-aspheric laser beam shapers,” Proc. SPIE **8490**, 849003 (2012). [CrossRef]

**14. **R. H. Estes and E. R. Lancaster, “Two-point Taylor series expansions,” NASA TMX-55759 (N67-23965) (1966).

**15. **Y. Cai, X. Lu, and Q. Lin, “Hollow Gaussian beams and their propagation properties,” Opt. Lett. **28**, 1084–1086 (2003). [CrossRef] [PubMed]

**16. **F. Duerr, Y. Meuret, and H. Thienpont, “Potential benefits of free-form optics in on-axis imaging applications with high aspect ratio,” Opt. Express **21**, 31072–31081 (2013). [CrossRef]

**17. **F. Duerr, Y. Meuret, and H. Thienpont, “Tailored free-form optics with movement to integrate tracking in concentrating photovoltaics,” Opt. Express **21**, A401–A411 (2013). [CrossRef] [PubMed]