Single freeform surface design for prescribed input wavefront and target   irradiance

Christoph B\"osel; Herbert Gross

arXiv:1705.04540·physics.optics·August 10, 2017

Single freeform surface design for prescribed input wavefront and target irradiance

Christoph B\"osel, Herbert Gross

PDF

TL;DR

This paper presents a novel PDE-based method for designing single freeform optical surfaces capable of transforming complex input wavefronts and irradiance distributions into desired target patterns, enhancing beam shaping flexibility.

Contribution

It introduces a PDE formulation for general zero-étendue sources and a numerical solution approach, enabling efficient design of freeform surfaces for complex beam shaping tasks.

Findings

01

Successfully designed freeform surfaces for various complex beam shaping scenarios

02

Demonstrated the method's effectiveness beyond paraxial approximation

03

Achieved accurate redistribution of diverse input wavefronts and irradiance patterns

Abstract

In beam shaping applications, the minimization of the number of necessary optical elements for the beam shaping process can benefit the compactness of the optical system and reduce its cost. The single freeform surface design for input wavefronts, which are neither planar nor spherical, is therefore of interest. In this work, the design of single freeform surfaces for a given zero-\'etendue source and complex target irradiances is investigated. Hence, not only collimated input beams or point sources are assumed. Instead, a predefined input ray direction vector field and irradiance distribution on a source plane, which has to be redistributed by a single freeform surface to give the predefined target irradiance, is considered. To solve this design problem, a partial differential equation (PDE) or PDE system, respectively, for the unknown surface and its corresponding ray mapping is…

Figures10

Click any figure to enlarge with its caption.

Tables3

Table 1. Table 1: Comparison of Δ I T Δ subscript 𝐼 𝑇 \Delta I_{T} for example “lake” with a collimated input beam.

	$z^{\infty} (x, y)$	$z (x, y)$
$r m s_{Δ I_{T}}$	$9.2288 \cdot 10^{- 6}$	$2.3472 \cdot 10^{- 6}$
$p v_{Δ I_{T}}$	$1.6158 \cdot 10^{- 4}$	$3.6349 \cdot 10^{- 5}$
$c o r r_{I_{T}}$	$0.5750$	$0.9607$
$η$	$85.66 %$	$99.98 %$

Table 2. Table 2: Comparison of Δ I T Δ subscript 𝐼 𝑇 \Delta I_{T} for example “lake” with a point source.

	$z^{\infty} (x, y)$	$z (x, y)$
$r m s_{Δ I_{T}}$	$4.1373 \cdot 10^{- 6}$	$2.1480 \cdot 10^{- 6}$
$p v_{Δ I_{T}}$	$9.0955 \cdot 10^{- 5}$	$3.1909 \cdot 10^{- 5}$
$c o r r_{I_{T}}$	$0.8889$	$0.9668$
$η$	$99.62 %$	$99.55 %$

Table 3. Table 3: Comparison of Δ I T Δ subscript 𝐼 𝑇 \Delta I_{T} for example “lake” with an astigmatic input wavefront

	$z^{\infty} (x, y)$	$z (x, y)$
$r m s_{Δ I_{T}}$	$3.2072 \cdot 10^{- 6}$	$2.7386 \cdot 10^{- 6}$
$p v_{Δ I_{T}}$	$4.2670 \cdot 10^{- 5}$	$3.9421 \cdot 10^{- 5}$
$c o r r_{I_{T}}$	$0.9252$	$0.9464$
$η$	$99.23 %$	$99.20 %$

Equations57

\hat{s}_{1} = (\hat{s}_{1})_{x} (\hat{s}_{1})_{y} (\hat{s}_{1})_{z},

\hat{s}_{1} = (\hat{s}_{1})_{x} (\hat{s}_{1})_{y} (\hat{s}_{1})_{z},

d e t (\nabla u (x, y)) I_{T} (u (x, y)) = I_{S} (x, y),

d e t (\nabla u (x, y)) I_{T} (u (x, y)) = I_{S} (x, y),

\hat{s}_{2} = n \hat{s}_{1} + {- n \cdot \hat{n} \cdot \hat{s}_{1} + 1 - n^{2} [1 - (\hat{n} \cdot \hat{s}_{1})^{2}]} \hat{n},

\hat{s}_{2} = n \hat{s}_{1} + {- n \cdot \hat{n} \cdot \hat{s}_{1} + 1 - n^{2} [1 - (\hat{n} \cdot \hat{s}_{1})^{2}]} \hat{n},

\hat{s}_{2} = \hat{s}_{1} - 2 (\hat{n} \cdot \hat{s}_{1}) \hat{n}

s_{1} = x_{S} - x y_{S} - y z (x_{S}, y_{S}) - z_{0}, s_{3} = u_{x} (x, y) - x u_{y} (x, y) - y z_{T} - z_{0},

s_{1} = x_{S} - x y_{S} - y z (x_{S}, y_{S}) - z_{0}, s_{3} = u_{x} (x, y) - x u_{y} (x, y) - y z_{T} - z_{0},

n = - \partial_{x_{S}} z (x_{S}, y_{S}) - \partial_{y_{S}} z (x_{S}, y_{S}) 1,

s_{1} + λ \cdot \hat{s}_{2} = s_{3} \Rightarrow λ = \frac{z _{T} - z ( x _{S} , y _{S} )}{( s ^ _{2} ) _{z}} .

s_{1} + λ \cdot \hat{s}_{2} = s_{3} \Rightarrow λ = \frac{z _{T} - z ( x _{S} , y _{S} )}{( s ^ _{2} ) _{z}} .

u_{x} (x, y) - x = λ \cdot (\hat{s}_{2})_{x} + (s_{1})_{x},

u_{x} (x, y) - x = λ \cdot (\hat{s}_{2})_{x} + (s_{1})_{x},

u_{y} (x, y) - y = λ \cdot (\hat{s}_{2})_{y} + (s_{1})_{y},

\frac{x _{S} - x}{z ( x _{S} , y _{S} ) - z _{0}} = \frac{( s ^ _{1} ) _{x}}{( s ^ _{1} ) _{z}} \Rightarrow x_{S} = \frac{( s ^ _{1} ) _{x}}{( s ^ _{1} ) _{z}} [z (x_{S}, y_{S}) - z_{0}] + x,

\frac{x _{S} - x}{z ( x _{S} , y _{S} ) - z _{0}} = \frac{( s ^ _{1} ) _{x}}{( s ^ _{1} ) _{z}} \Rightarrow x_{S} = \frac{( s ^ _{1} ) _{x}}{( s ^ _{1} ) _{z}} [z (x_{S}, y_{S}) - z_{0}] + x,

\frac{y _{S} - y}{z ( x _{S} , y _{S} ) - z _{0}} = \frac{( s ^ _{1} ) _{y}}{( s ^ _{1} ) _{z}} \Rightarrow y_{S} = \frac{( s ^ _{1} ) _{y}}{( s ^ _{1} ) _{z}} [z (x_{S}, y_{S}) - z_{0}] + y,

\partial_{x_{S}} z (x_{S}, y_{S}) = \frac{( \partial _{y} y _{S} ) [ \partial _{x} z ( x _{S} , y _{S} ) ] - ( \partial _{x} y _{S} ) [ \partial _{y} z ( x _{S} , y _{S} )]}{( \partial _{x} x _{S} ) \cdot ( \partial _{y} y _{S} ) - ( \partial _{x} y _{S} ) \cdot ( \partial _{y} x _{S} )},

\partial_{x_{S}} z (x_{S}, y_{S}) = \frac{( \partial _{y} y _{S} ) [ \partial _{x} z ( x _{S} , y _{S} ) ] - ( \partial _{x} y _{S} ) [ \partial _{y} z ( x _{S} , y _{S} )]}{( \partial _{x} x _{S} ) \cdot ( \partial _{y} y _{S} ) - ( \partial _{x} y _{S} ) \cdot ( \partial _{y} x _{S} )},

\partial_{y_{S}} z (x_{S}, y_{S}) = \frac{- ( \partial _{y} x _{S} ) [ \partial _{x} z ( x _{S} , y _{S} )] + ( \partial _{x} x _{S} ) [ \partial _{y} z ( x _{S} , y _{S} )]}{( \partial _{x} x _{S} ) \cdot ( \partial _{y} y _{S} ) - ( \partial _{x} y _{S} ) \cdot ( \partial _{y} x _{S} )} .

z_{S} (x, y) \equiv z (x_{S}, y_{S}),

z_{S} (x, y) \equiv z (x_{S}, y_{S}),

[(\partial_{x} u_{x}) (\partial_{y} u_{y}) - (\partial_{x} u_{y}) (\partial_{y} u_{x})] I_{T} (u_{x}, u_{y}) = I_{S} (x, y),

[(\partial_{x} u_{x}) (\partial_{y} u_{y}) - (\partial_{x} u_{y}) (\partial_{y} u_{x})] I_{T} (u_{x}, u_{y}) = I_{S} (x, y),

u_{x} (x, y) - x = [z_{T} - z_{S} (x, y)] \frac{( s ^ _{2} ) _{x}}{( s ^ _{2} ) _{z}} + [z_{S} (x, y) - z_{0}] \frac{( s ^ _{1} ) _{x}}{( s ^ _{1} ) _{z}},

u_{y} (x, y) - y = [z_{T} - z_{S} (x, y)] \frac{( s ^ _{2} ) _{y}}{( s ^ _{2} ) _{z}} + [z_{S} (x, y) - z_{0}] \frac{( s ^ _{1} ) _{y}}{( s ^ _{1} ) _{z}},

A \partial_{xx} z_{S} + 2 B \partial_{x y} z_{S} + C \partial_{y y} z_{S} + [(\partial_{xx} z_{S}) (\partial_{y y} z_{S}) - (\partial_{x y} z_{S})^{2}] = D .

A \partial_{xx} z_{S} + 2 B \partial_{x y} z_{S} + C \partial_{y y} z_{S} + [(\partial_{xx} z_{S}) (\partial_{y y} z_{S}) - (\partial_{x y} z_{S})^{2}] = D .

d (I_{S}, I_{T})^{2} = u \in M in f \int ∣ u (x) - x ∣^{2} I_{S} (x) d x,

d (I_{S}, I_{T})^{2} = u \in M in f \int ∣ u (x) - x ∣^{2} I_{S} (x) d x,

\frac{n ( x _{S} , y _{S} )}{( n ( x _{S} , y _{S} ) ) _{z}} = \frac{n _{1} s ^ _{1} - n _{2} s ^ _{2}}{( n _{1} s ^ _{1} - n _{2} s ^ _{2} ) _{z}} = \nabla_{S} (z - z_{S} (x, y))

\frac{n ( x _{S} , y _{S} )}{( n ( x _{S} , y _{S} ) ) _{z}} = \frac{n _{1} s ^ _{1} - n _{2} s ^ _{2}}{( n _{1} s ^ _{1} - n _{2} s ^ _{2} ) _{z}} = \nabla_{S} (z - z_{S} (x, y))

\partial_{x_{S}} z_{S}^{\infty} (x, y) = - \frac{n _{1} ( s ^ _{1} ) _{x} - n _{2} ( s ^ _{2} ) _{x}}{n _{1} ( s ^ _{1} ) _{z} - n _{2} ( s ^ _{2} ) _{z}} \equiv n_{x} (x, y),

\partial_{x_{S}} z_{S}^{\infty} (x, y) = - \frac{n _{1} ( s ^ _{1} ) _{x} - n _{2} ( s ^ _{2} ) _{x}}{n _{1} ( s ^ _{1} ) _{z} - n _{2} ( s ^ _{2} ) _{z}} \equiv n_{x} (x, y),

\partial_{y_{S}} z_{S}^{\infty} (x, y) = - \frac{n _{1} ( s ^ _{1} ) _{y} - n _{2} ( s ^ _{2} ) _{y}}{n _{1} ( s ^ _{1} ) _{z} - n _{2} ( s ^ _{2} ) _{z}} \equiv n_{y} (x, y) .

\partial_{x} z_{S}^{\infty} (x, y) =

\partial_{x} z_{S}^{\infty} (x, y) =

\frac{n _{x} ( x , y ) + [ n _{x} ( x , y ) \partial _{x} ( \frac{( s ^ _{1} ) _{x}}{( s ^ _{1} ) _{z}} ) + n _{y} ( x , y ) \partial _{x} ( \frac{( s ^ _{1} ) _{y}}{( s ^ _{1} ) _{z}} ) ] ( z _{S}^{\infty} ( x , y ) - z _{0} )}{1 - n _{x} ( x , y ) \frac{( s ^ _{1} ) _{x}}{( s ^ _{1} ) _{z}} - n _{y} ( x , y ) \frac{( s ^ _{1} ) _{y}}{( s ^ _{1} ) _{z}}},

\partial_{y} z_{S}^{\infty} (x, y) =

\frac{n _{y} ( x , y ) + [ n _{x} ( x , y ) \partial _{y} ( \frac{( s ^ _{1} ) _{x}}{( s ^ _{1} ) _{z}} ) + n _{y} ( x , y ) \partial _{y} ( \frac{( s ^ _{1} ) _{y}}{( s ^ _{1} ) _{z}} ) ] ( z _{S}^{\infty} ( x , y ) - z _{0} )}{1 - n _{x} ( x , y ) \frac{( s ^ _{1} ) _{x}}{( s ^ _{1} ) _{z}} - n _{y} ( x , y ) \frac{( s ^ _{1} ) _{y}}{( s ^ _{1} ) _{z}}}

u (x, y)

u (x, y)

z_{S} (x, y)

\partial_{x} (Δ z_{S}) \to \frac{1}{2Δ x} [(Δ z_{S})_{i; j + 1} - (Δ z_{S})_{i; j - 1}],

\partial_{x} (Δ z_{S}) \to \frac{1}{2Δ x} [(Δ z_{S})_{i; j + 1} - (Δ z_{S})_{i; j - 1}],

\partial_{y} (Δ z_{S}) \to \frac{1}{2Δ y} [(Δ z_{S})_{i + 1; j} - (Δ z_{S})_{i - 1; j}]

\partial_{x} (Δ z_{S}) \to - \frac{1}{2Δ x} [3 (Δ z_{S})_{i; j + 2} - 4 (Δ z_{S})_{i; j + 1} + (Δ z_{S})_{i; j}],

\partial_{x} (Δ z_{S}) \to - \frac{1}{2Δ x} [3 (Δ z_{S})_{i; j + 2} - 4 (Δ z_{S})_{i; j + 1} + (Δ z_{S})_{i; j}],

\partial_{y} (Δ z_{S}) \to - \frac{1}{2Δ y} [3 (Δ z_{S})_{i + 2; j} - 4 (Δ z_{S})_{i + 1; j} + (Δ z_{S})_{i; j}],

Δ u_{x} (- 0.5, y) = Δ u_{x} (0.5, y) = 0, y \in [- 0.5, 0.5],

Δ u_{x} (- 0.5, y) = Δ u_{x} (0.5, y) = 0, y \in [- 0.5, 0.5],

Δ u_{y} (x, - 0.5) = Δ u_{y} (x, 0.5) = 0, x \in [- 0.5, 0.5] .

\hat{s}_{1} = 001

\hat{s}_{1} = 001

\hat{s}_{1} = \frac{1}{x ^{2} + y ^{2} + z _{0}^{2}} x y z_{0} .

\hat{s}_{1} = \frac{1}{x ^{2} + y ^{2} + z _{0}^{2}} x y z_{0} .

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

Single freeform surface design for prescribed input wavefront and target irradiance

Christoph Bösel

[email protected]

Friedrich-Schiller-Universität Jena, Institute of Applied Physics, Albert-Einstein-Str. 15, 07745 Jena, Germany

Herbert Gross

Friedrich-Schiller-Universität Jena, Institute of Applied Physics, Albert-Einstein-Str. 15, 07745 Jena, Germany

Fraunhofer-Institut für Angewandte Optik und Feinmechanik, Albert-Einstein-Str. 7, 07745 Jena, Germany

Abstract

In beam shaping applications, the minimization of the number of necessary optical elements for the beam shaping process can benefit the compactness of the optical system and reduce its cost. The single freeform surface design for input wavefronts, which are neither planar nor spherical, is therefore of interest. In this work, the design of single freeform surfaces for a given zero-étendue source and complex target irradiances is investigated. Hence, not only collimated input beams or point sources are assumed. Instead, a predefined input ray direction vector field and irradiance distribution on a source plane, which has to be redistributed by a single freeform surface to give the predefined target irradiance, is considered. To solve this design problem, a partial differential equation (PDE) or PDE system, respectively, for the unknown surface and its corresponding ray mapping is derived from energy conservation and the ray-tracing equations. In contrast to former PDE formulations of the single freeform design problem, the derived PDE of Monge-Ampère type is formulated for general zero-étendue sources in cartesian coordinates. The PDE system is discretized with finite differences and the resulting nonlinear equation system solved by a root-finding algorithm. The basis of the efficient solution of the PDE system builds the introduction of an initial iterate constuction approach for a given input direction vector field, which uses optimal mass transport with a quadratic cost function. After a detailed description of the numerical algorithm, the efficiency of the design method is demonstrated by applying it to several design examples. This includes the redistribution of a collimated input beam beyond the paraxial approximation, the shaping of point source radiation and the shaping of an astigmatic input wavefront into a complex target irradiance distribution.

PACS numbers

42.15.-i, 42.15.Eq.

I Introduction

The development of numerical design algorithms for the calculation of freeform surfaces for irradiance control in the geometrical optics limit is of great interest in current research. The task of these methods is thereby the calculation of freeform surfaces, which redistribute a given light source into a required target irradiance pattern. In particular the calculation of continuous freeform surfaces, which is considerd here, is of importance for manufacturability and the reduction of possible diffraction effects. Especially the irradiance control with continuous single freeform surfaces for input beams with planar or nonspherical wavefronts lacks of detailed discussion in the literature. Such design algorithms can potentially reduce the number of necessary optical elements or increase the efficiency in practical applications and are investigated in this work.

In principle, the freeform design for shaping of ideal sources can be divided into three major groups. The first group consists of design approaches using Oliker’s supporting quadrics method (SQM) Oliker02_1 ; Fou10_1 ; Mic11_1 ; Oliker17_1 , which is arguably the theoretically best justified, most reliable freeform design method and usable in complex optical systems Mic11_1 . As it was shown by Oliker Oliker02_1 , the geometric properties of quadric surfaces can be used effectively for the freeform design process. The shaping of e.g. a point source with a freeform mirror can be done by placing the point source in the common focal point of a number of ellipsoids. By changing the focal parameters of each ellipsoid appropriately, a freeform mirror, which gives the required irradiance distribution on the target area, can be designed by a unification of the ellipsoids. Due to the high computational effort, necessary for high resolution target irradiances and large number of quadrics, the major difficulty, considering the SQM, is the development of efficient algorithms.

The second group of design approaches consists of methods, which are modelling the freeform design problem by a highly nonlinear partial differential equation (PDE) of Monge-Ampère type for point sources ries02 ; Wu13_1 ; Brix15_1 and collimated beams Wu13_2 ; Prins14_1 . The main problem is thereby the development of numerical algorithms to solve the PDE. An advantage of the PDE methods over the SQM is the lower computational effort, enabling them to handle high resolution target irradiances. The major drawback of the PDE methods is the lack of theoretical justification and reliability as pointed out by Oliker Oliker17_1 .

The third group are ray-mapping or source-target mapping methods, respectively Oli13_1 ; Brun13_1 ; Feng13_2 ; Ma15_1 ; Boe16_1 ; Gan17_1 . These are based on the calculation of a ray mapping between the source and the target irradiance distribution an the subsequent construction of the freeform surface from the mapping. The major difficulty of these methods is the finding of an integrable mapping, enabling the design of a freeform surface, which is on the one hand continuous and on the other hand maps the source onto the target distribution. Some of the ray-mapping methods Brun13_1 ; Feng13_2 ; Boe16_1 explore a mapping calculated from optimal mass transport (OMT) theory with a quadratic cost function. This special mapping solves the collimated beam shaping problem in a paraxial approximation Boe16_1 , but is in general not integrable. As we will demonstrate in this work, the quadratic cost function mapping is nevertheless of great usefulness for the freeform design process for more general optical configurations.

In this work, we want to generalize the design approach for collimated beam shaping with continuous single freeform surfaces presented in Ref. Boe16_1 to general input wavefronts of ideal sources, whereby, contrary to former PDE formulations, parallel input beams or point sources are not assumed. Instead, we derive a system of PDE’s and a Monge-Ampère equation, respectively, for a predefined input direction vector field and present a detailed numerical algorithm to solve the corresponding problem. As indicated in the previous paragraph, the basis of the numerical algorithm builds thereby a ray mapping, derived from OMT with a quadratic cost function, which will serve as an initial iterate for solving the derived PDE system. As we will demonstrate explicitely, the design algorithm extends the work of Ref. Boe16_1 firstly to collimated beam shaping in the nonparaxial regime, secondly to the irradiance control of point sources, and thirdly, to other input wavefront shapes.

To introduce the design algorithm, this work is structured as follows. In section II and III, a PDE system is derived from energy conservation and the ray-tracing equations for single freeform lens and mirror surfaces. For the numerical solution of the PDE system, a method for the calculation of an appropriate initial iterate is introduced in section IV. In section V, the numerical algorithm for determining the desired freeform surface is presented and in section VI several design examples, illustrating the effectiveness and usefulness of the design algorithm, are discussed.

II Geometry

In the following, the problem of mapping an incoming beam, defined by the ray direction vector field

[TABLE]

and the irradiance distribution $I_{S}(x,y)$ on the plane $z=z_{0}$ , onto a predefined output irradiance distribution $I_{T}(x,y)$ on the target plane $z=z_{T}$ is considered (see Fig. 1).

The goal is thereby the calculation of a continuous reflecting or refracting surface, which enables the mapping of the source onto the target distribution. In a geometrical optics framework, this problem is fully governed by two equations. This is on the one hand the Jacobian equation

[TABLE]

describing the relation of the input and output irradiance by the ray mapping or coordinate transformation $\mathbf{u}(x,y)$ , respectively, and on the other hand the ray tracing equations for refractive and reflective surfaces

[TABLE]

with the surface unit normal vector field $\mathbf{\hat{n}}$ . Equation (3) relates the input coordinates $(x,y)$ and the target coordinates $\mathbf{u}(x,y)$ by the unknown surface $z(x,y)$ . Unit vector fields are hereby denoted by the “hat”. To determine the unknown surface, we want to derive a PDE system for the mapping $\mathbf{u}(x,y)$ and the surface $z(x,y)$ . This can be done by expressing the considered geometry, according to Fig. 1, through the vector fields

[TABLE]

and by calculating the intersection points $\mathbf{u}(x,y)$ with the target plane $z=z_{T}$ from the relation

[TABLE]

Hereby $x_{S},y_{S}$ and $z(x_{S},y_{S})$ are the intersection coordinates of the input vector field $\mathbf{s_{1}}$ with the surface $z(x,y)$ . From Eq. (5) we then directly get

[TABLE]

which define a relation between the surface $z(x,y)$ , the surface gradient and the target coordinates. For the case of collimated input beams, analogues equations were derived by Wu et al. Wu13_2 and Oliker Oli14_1 . Together with Eq. (2), Eq. (6) builds a PDE system for the unknown functions $u_{x}(x,y),u_{y}(x,y)$ and $z(x,y)$ . According to Eq. (4), the remaining difficulty is thereby the dependency of Eq. (6) on the surface intersection points $(x_{S},y_{S})$ , which are determined by the surface $z(x,y)$ and therefore unknown. A similar problem arising in their work on the reconstruction of an optical surface from experimental ray data Rub01_1 , was solved by Rubinstein and Wolansky using the geometrical relation

[TABLE]

between the given ray directions defined in Eq. (1) and the vector field $\mathbf{s}_{1}$ in Eq. (4). By applying this coordinate transformation to the coordinates $x_{S}$ , $y_{S}$ and the gradient $\nabla_{S}z(x_{S},y_{S})=(\partial_{x_{S}}z(x_{S},y_{S}),\partial_{y_{S}}z(x_{S},y_{S}))$ occuring in Eq. (4), Eq. (6) can be expressed by the initial coordinates $(x,y)$ on the plane $z=z_{0}$ . After using the chain rule it follows

[TABLE]

The remaining derivatives of $x_{S}$ and $y_{S}$ can be evaluated with Eq. (7). With these transformations and the definition

[TABLE]

Eq. (6) is fully expressed in terms of the $x$ - $y$ -coordinates. Hence, the Jacobian equation (2) and Eq. (6) build a system of PDEs for the unknown functions $\mathbf{u}(x,y)$ and $z_{S}(x,y)$ . After solving the equation system for $z_{S}(x,y)$ , the coordinates $(x_{S},y_{S})$ can than be calculated from Eq. (7), which gives the surface $z(x,y)$ on the scattered $x_{S}$ - $y_{S}$ grid.

III PDE system and Monge-Ampère equation

In the previous section it was shown that the single freeform design problem for zero-étendue sources is governed by the PDE system

[TABLE]

whereby $\mathbf{\hat{s}_{1}}$ is the predefined input ray direction vector field and $\mathbf{\hat{s}_{2}}$ is defined according to section II in terms of $z_{S}(x,y)$ , $\partial_{x}z_{S}(x,y)$ and $\partial_{y}z_{S}(x,y)$ .

For a given source and target distribution $I_{S}(x,y)$ and $I_{T}(x,y)$ , and a given input ray direction vector field $\mathbf{\hat{s}_{1}}$ the equation system in Eq. (10) has to be solved simultaneously for the surface $z_{S}(x,y)$ and the ray mapping $\mathbf{u}(x,y)$ .

Equation (10) can be further simplified by plugging $u_{x}(x,y)$ and $u_{y}(x,y)$ into the Jacobian equation to obtain a nonlinear second order PDE of Monge-Ampère type for the function $z_{S}(x,y)$ :

[TABLE]

Hereby, $A,B,C,D$ are functions of $x,y,z_{S}(x,y),\partial_{x}z_{S}(x,y)$ and $\partial_{y}z_{S}(x,y)$ , which can be calculated straightforeward.

Compared to the Eq. system (10), the Eq. (11) has numerically the advantage of a lesser number of design variables. Despite this advantage, the numerical method presented in the following will nevertheless be based on the description in Eq. (10). In contrast to the Monge-Ampère equation in Eq. (11), we will therefore omit the discretization of second order derivatives. In addition to that, the boundary conditions can be directly controlled through the mapping $\mathbf{u}(x,y)$ by requiring that boundary points of the support of $I_{S}(x,y)$ are mapped onto boundary points of the support of $I_{T}(x,y)$ . This is contrary to the Monge-Ampère equation for which the same goal has to be accomplished through the surface gradients $\nabla z_{S}(x,y)$ at the boundary of the surface, which leads to the dependency of the boundary conditions on the input field $\mathbf{\hat{s}}_{1}$ . Furthermore, the properties of the Jacobian equation and the mapping equations in Eq. (10) can be controlled separately.

IV Initial iterate

The construction of freeform surfaces for complex target irradiances from nonlinear PDEs, like Eqs. (10) or (11) is usually done by discretizatizing the PDEs and solving the resulting nonlinear equation system by an iterative process. A high resolution target irradiance distribution and therefore large number of unknowns requires in general an appropriate initial iterate to ensure (fast) convergence. Especially for discretization schemes, which are not provably convergent, those initial iterates are of great importance. For the published PDE methods for single freeform surface design an inital iterate producing a uniform illumation on the target plane or other simple surface shapes are usually chosen.

Here, we present a different method to construct an initial iterate for a given input direction vector field $\mathbf{\hat{s}_{1}}$ . It is based on the OMT with the quadratic cost function and inspired by Ref. Rub01_1 . Following the argumentation from Ref. Rub01_1 , an unknown surface can be experimentally reconstructed from the known input ray directions and measured target plane intersection points of the rays send out from $(x,y)$ . Since the target points are not given for the freeform design problem considered here, the approach from Ref. Rub01_1 can not be directly applied. Nevertheless, as we will demonstrate in the following, the approach is of great usefullness to generate an initial iterate for the freeform surface design.

IV.1 Initial mapping

For the single freeform surface design process a frequently used mapping $\mathbf{u}(x,y)$ is calculated from OMT by considering Eq. (2) together with the quadratic cost function condition

[TABLE]

which we will refer to in the following by $\mathbf{u^{\infty}}(x,y)$ . $M$ in Eq. (12) denotes the set of energy preserving mappings. This is arguably the most studied OMT problem in numerical mathematics and provides us therefore with efficient algorithms for the mapping calculation. As it was demonstrated in several publications, the corresponding mapping $\mathbf{u^{\infty}}(x,y)$ in general will not generate a continuous surface $z(x,y)$ fulfilling the required design goal of mapping the given input $I_{S}(x,y)$ onto the desired output irradiance $I_{T}(x,y)$ . However, due to its properties, the mapping $\mathbf{u^{\infty}}(x,y)$ appears to be a good candidate for creating an appropriate initial iterate. According to Eq. (12), the quadratic cost function can in principle be thought of as describing a ray propagation in a homogeneous medium Rub04_1 . Since a continuous surface can be interpreted as a perturbation of this free propagation, it can be expected that the mapping $\mathbf{u^{\infty}}(x,y)$ serves as a reasonable initial iterate for solving Eq. (10). From a numerical point of view it is also helpful that $\mathbf{u^{\infty}}(x,y)$ already fulfills the Jacobian equation in Eq. (10), which therefore basically constrains all numerical corrections to $\mathbf{u^{\infty}}(x,y)$ to the space of energy preserving mappings. Hence, we will use the mapping defined by Eqs. (2) and (12) as an initial mapping for solving Eq. (10). As in the previous publications Boe16_1 ; Boe17_1 , we will make use of the algorithm by Sulman et al. Sul11_1 for the examples in section VI for demonstrational purposes.

IV.2 Initial surface

After its calculation, the mapping $\mathbf{u}^{\infty}(x,y)$ can be used to construct the initial surface $z_{S}^{\infty}(x,y)$ as demonstrated in the following. The resulting initial surface construction equations are thereby a generalization of the surface construction equations used in Ref. Boe16_1 to noncollimated beams.

The derivation can be done by using the law of refraction/reflection

[TABLE]

and applying again the coordinate transformation defined in Eq. (7) to invert the gradient by using

[TABLE]

After using the chain rule, this gives the differential equations Rub01_1

[TABLE]

for the initial surface $z_{S}^{\infty}(x,y)$ . For a given input direction vector field, defined by Eq. (1), and given target points $\mathbf{u}^{\infty}(x,y)$ , Eq. (15) can be integrated along an arbitrary path on the support of $I_{S}(x,y)$ to generate the initial surface $z_{S}^{\infty}(x,y)$ .

Together with the initial mapping $\mathbf{u}^{\infty}(x,y)$ , the surface $z_{S}^{\infty}(x,y)$ builds an initial iterate for solving the PDE system (10).

We want to remark that the initial surface construction approach, presented in this section can also be used for the Monge-Ampère equation in Eq. (11). To implement the boundary conditions in that case, the surface boundary normal vectors defined by Eq. (15) have to be chosen in a way that the boundary of $I_{S}(x,y)$ is mapped onto the boundary of $I_{T}(x,y)$ .

V Numerics

V.1 Discretization

After the calculation of the initial iterate, the PDE system (10) has to be solved. Hence, we have to discretize (10) and solve the resulting nonlinear equation system by a root-finding algorithm. To do that, we follow Ref. Boe17_1 by introducing the corrections $\Delta\mathbf{u}(x,y)$ and $\Delta z_{S}(x,y)$ to $\mathbf{u}^{\infty}(x,y)$ and $z^{\infty}(x,y)$ :

[TABLE]

and redefine the initial iterates for the solving process to be $\Delta\mathbf{u}(x,y)=0$ and $\Delta z(x,y)=0$ . Accordingly, Eq. (10) becomes a PDE system for $\Delta\mathbf{u}(x,y)$ and $\Delta z_{S}(x,y)$ . To solve the PDE system, we discretize the functions $\Delta u_{x}(x,y),\Delta u_{y}(x,y),\Delta z_{S}(x,y)$ and their derivatives in Eq. (10) by central finite differences

[TABLE]

for the inner points $i;j=2,...,N-1$ and second order finite differences

[TABLE]

for the boundary points. Additionally, we apply boundary conditions by requiring that the boundary points of the source distribution $I_{S}(x,y)$ are mapped onto the boundary points of the target distribution $I_{T}(x,y)$ . For the mapping of square supports, which is demonstrated in section 6, the boundary condition reads therefore:

[TABLE]

After the discretization of Eq. (10), the resulting nonlinear Eq. system can then be solved by a root-find algorithm for the $3\cdot N^{2}$ unknown values $(\Delta u_{x})_{i;j},(\Delta u_{y})_{i;j}$ and $(\Delta z_{S})_{i;j}$ . The solution will lead to a discretized surface $z_{S}(x,y)$ from which, according to Eq. (7), the scattered grid points $(x_{S},y_{S})$ can be calculated. As in Ref. Boe17_1 , we use MATLAB’s fsolve() function with the trust-region-reflective solver for the solution of the nonlinear Eq. system (10).

V.2 Design algorithm

The design algorithm developed in the previous sections can be summarized as follows:

Define the input ray direction vector field $\mathbf{\hat{s}_{1}}$ and the input irradiance $I_{S}(x,y)$ on $z=z_{0}$ and the output irradiance $I_{T}(x,y)$ on $z=z_{T}$ . 2. 2.

Calculate the initial map $\mathbf{u}^{\infty}(x,y)$ between $I_{S}(x,y)$ and $I_{T}(x,y)$ . 3. 3.

Integrate Eq. (15) along an arbitrary path on the support of $I_{S}(x,y)$ to get the initial surface $z_{S}^{\infty}(x,y)$ . 4. 4.

Discretize the PDE system in Eq. (10) to get a nonlinear equation system. 5. 5.

Apply the boundary conditions and solve the nonlinear equation system with $\mathbf{u}^{\infty}(x,y)$ and $z_{S}^{\infty}(x,y)$ as the initial iterate to get $\mathbf{u}(x,y)$ and $z_{S}(x,y)$ .

As mentioned above, for a general input vector field $\mathbf{\hat{s}_{1}}$ , the presented design algorithm gives us the freeform surface $z(x,y)$ on the scattered grid points $(x_{S},y_{S})$ . The quality of the surface $z_{S}(x,y)$ is thereby already characterized by Eq. (10), which is equivalent to a ray-tracing calculation. For manufacturing processes it might still be useful to evaluate the surface quality in a ray-tracing simulation. We therefore interpolate the scattered data points $z(x_{S},y_{S})$ onto a regular grid and import it to a ray-tracing software. The interpolation can be done by scattered data interpolation with e.g. MATLAB 2015b’s griddata() function. However, in our experience this function has two major drawbacks for the interpolation of freeform surfaces. This is its inability (for large data sets) to extrapolate outside of the convex hull of the grid points, which leads to inaccuracies at the boundary of the freeform surface in form of zig-zag patterns at the target irradiance boundary. Furthermore, the interpolation with griddata() leads to oscillations in the freeform surface, depending on the stepsize of the interpolation grid and the interpolation method.

To overcome this limitation, we therefore use the RegularizeData3D() function from the MATLAB file exchange website MatFE_1 . This function provides the possibility to choose how smooth the interpolated data is fitted to the scattered data. In our examples, we made good experience with a small smootheness factor for low noise data. Additionally, RegularizeData3D() extrapolates outside the convex hull of the scattered data grid, which may also lead to slight inaccuracies at the boundary of the freeform surface, which lower the energy efficiency, as we will point out in the example section VI.

Since the initial mapping $\mathbf{u}^{\infty}(x,y)$ defined by Eqs. (2) and (12) is in general not integrable, the integration of Eq. (15) will be path dependent and therefore lead to different initial surfaces $z_{S}^{\infty}(x,y)$ for different integration paths. Nevertheless, in our experience the path dependency has no significant influence on the suitability for using $z_{S}^{\infty}(x,y)$ as an initial iterate.

VI Examples

In this section, we want to demonstrate basically three things. Firstly, the design algorithm presented in section V extends the work of Ref. Boe16_1 on collimated beam shaping with single freeform surfaces to the nonparaxial regime. Secondly, the design algorithm can be applied to the design of single freeform surfaces for shaping of point light sources, which is frequently discussed in literature. And thirdly, we want to demonstrate the applicability to a nonstandard design example by dealing with an input wavefront, which is neither planar nor spherical. This is done in the following three subsections. All the examples in this sections were calculated on an Intel Core i3 at $2\times 2.4$ Ghz with $16$ GB RAM. To make the ray-tracing results comparable, we used for each of the examples the target distribution “lake” (see Fig. 3(b)) with a resolution of $250\times 250$ pixels. For the calculation of the initial surface by Eq. (15), we used MATLAB’s ode45 solver with tolerances of $10^{-10}$ , which is in our experience necessary for a legitimate comparison of the initial and final surface in a ray-tracing simulation. For the root finding of Eq. (10) by fsolve() tolerances of $10^{-3}$ were used. For the examples with lenses refractive indices of $n=1.5$ were used. The ray-tracing simulations for the following examples are done with $200\cdot 10^{6}$ rays and detectors with $250\times 250$ pixels by using a self-programmed MATLAB ray-tracing toolbox. Since the design examples are calculated in a geometrical optics framework, we will use arbitrary units as length units and omit the explicit declaration of lengths by “a.u.”.

VI.1 Collimated beam shaping beyond the paraxial approximation

In Ref. Boe16_1 it was shown that the single freeform surface design problem for collimated input beams

[TABLE]

can be solved by the quadratic cost function ray mapping, defined by Eqs. (2) and (12), in the paraxial regime. Here, we demonstrate that the presented algorithm from section V directly extends the results from Ref. Boe16_1 to nonparaxial calculations. Since the quadratic cost function mapping solves the single freeform surface design problem exactly for infinite distances between the freeform surface and the target planeBoe16_1 , it can be argued analogously to Ref. Boe17_1 that the numerical convergence of the algorithm can theoretically be ensured by a stepwise reduction of the target plane distance $z_{T}$ . However, in our experience, this stepwise reduction is in a practical calculation not necessary to ensure numerical convergence.

From an application point of view, the extension to the nonparaxial regime is especially important for optical configurations in which the distance between the target plane and the freeform surface is small, the target area is strongly shifted or scaled, or if there are caustics between the target plane and the freeform surface. In the following, we therefore demonstrate the utility of the presented algorithm by designing a freeform mirror, which shifts and scales the target area compared to the position and size of the input beam. The corresponding system layout, considered in this section, can be seen in Fig. 2, showing a target area with a side length, which is $3$ times as large as the input beams side length of $1$ and shifted by $5$ . The distance between the freeform mirror and the target plane $z_{T}$ is approximately $5$ .

For the irradiance distribution, we choose the Gaussian distribution, which is represented in Fig. 3 (a) and the test image “lake” as the required target distribution (Fig. 3 (b)). According to section V, we first calculate the initial mapping between the input and output irradiances ( $320sec$ ), then the initial surface from Eq. (15) ( $68sec$ ) and finally solve the PDE system in Eq. (10) ( $208sec$ ). After the construction of the initial and final freeform mirror, we evaluate their quality by a ray-tracing simulation. The corresponding results and comparisons with the algorithm presented in Ref. Boe16_1 , which is equal to the initial iterate of the root finding, are shown in Fig. 4. Additionally, the $rms$ and the peak-to-valley of $\Delta I_{T}(x,y)=I_{T}(x,y)-I_{T,RT}(x,y)$ , and the correlation coefficient $corr_{I_{T}}$ between the predefined irradiance $I_{T}(x,y)$ and the irradiance from the ray tracing $I_{T,RT}(x,y)$ , and the energy efficiency $\eta$ are given in Table 1 . The results thereby show clearly the improvement of the presented design algorithm compared to Ref. Boe16_1 and the mapping from the quadratic cost function (12), respectively.

VI.2 Point source

Besides the collimated beam shaping with a single freeform surface, we want to apply the design algorithm to the other standard configuration considered in literature, which is the mapping of a point source onto a required irradiance distribution on a target plane $z_{T}$ . In the case of single freeform lens design, the point source is assumed to be in the medium or the rays emitted by the point source first hit a spherical surface, respectively. The ray directions are therefore first unaltered and then redirected by the freeform surface to give the predefined target distribution $I_{T}(x,y)$ .

Since we want to apply the algorithm from section V, we first have to project the point source onto the input plane $z=z_{0}$ , giving the distribution $I_{S}(x,y)$ . Subsequently, the mapping $\mathbf{u}^{\infty}(x,y)$ between the projected distribution $I_{S}(x,y)$ on $z=z_{0}$ and the target distribution $I_{T}(x,y)$ on $z=z_{T}$ can be calculated. For the point source, we choose a lambertian distribution with a maximum opening angle of $30$ degrees, which refers to the angle between the $z$ -axis and the corner points of the quadratic aperture and leads to the irradiance distribution $I_{S}(x,y)$ presented in Fig. 5(a).

The input vector field $\mathbf{\hat{s}_{1}}$ on $z=z_{0}$ is given by

[TABLE]

The placement of the plane $z=z_{0}$ relative to the point source does not matter. Depending on the value of $z_{0}$ only the mapping $\mathbf{u}^{\infty}(x,y)$ has to be scaled by a constant factor to create the predefined target distribution area size, which we choose to have a side length of $5$ . For the distance between the point source and the freeform surface we choose $z(0,0)=1$ , which enters the design process as an integration constant in Eq. (15), and for the distance between the point source and the target plane, we choose $z_{T}=10$ . The corresponding system layout can be seen in Fig. 6.

The calculation of the initial mapping $\mathbf{u}^{\infty}(x,y)$ took approximately $378sec$ , the initial surface $z_{S}^{\infty}(x,y)$ $183sec$ and the final surface $z_{S}(x,y)$ about $213sec$ . After the calculation of the freeform surface $z_{S}(x,y)$ by the presented design algorithm an additional complication, compared to the collimated beam shaping, arises since the freeform surface is calculated on the nonuniform $x_{S}$ - $y_{S}$ -grid. Fig. 5(b) shows every 4th point of the scattered grid points. The outermost points define the aperture of the freeform lens.

As described in section V, the function $z_{S}(x,y)=z(x_{S},y_{S})$ is interpolated and extrapolated onto a square, uniform grid by RegularizeData3D() to give $z(x,y)$ . For a fine boundary interpolation, we choose an interpolation grid of $500\times 500$ points, whereby the minimum and maximum $x$ and $y$ values correspond to the minimum and maximum $x_{S}$ and $y_{S}$ values (see Fig. 5(b)). The function values are then imported to a ray-tracing software. In Fig. 7 and Table 2 the results of the ray tracing are shown. From Fig. 7(a) and 7(c) it can be seen that the quadratic cost function in Eq. (12) provides a reasonable initial iterate for the root-finding process, as argued in section IV. Figure 7(b) and 7(c) and Table 2 show results, which are comparable to the example from the previous subsection , except for the slightly lowered efficiency $\eta$ due to the interpolation of $z_{S}(x,y)$ at the surface boundary.

From the viewpoint of manufacturing, it is usually preferred to have a plane instead of a spherical surface as the entrance surface of the lens. In the following section it is demonstrated that the design algorithm can be applied directly to these configurations.

VI.3 Astigmatic wavefront

In this section, we want to explore an example with a nonstandard input vector field $\mathbf{\hat{s}}_{1}$ . We therefore predefine a surface $z_{pre}(x,y)=\frac{1}{9}y^{2}-\frac{2}{3}x^{2}$ , deforming a collimated beam with a gaussian distribution (Fig. 9(a)) into a beam with an astigmatic wavefront. The corresponding system layout can be seen in Fig. 8.

As in the previous examples, we need the input vector field $\mathbf{\hat{s}}_{1}$ and the irradiance distribution on an intermediate plane $z=z_{0}$ . This plane is chosen to be $z_{0}=5$ and therefore after the refraction of the astigmatic beam by the plane entrance surface ( $z=3$ in Fig. 8(a)), and the aperture at $z=z_{0}$ is a square with a side length of $1$ . The vector field $\mathbf{\hat{s}}_{1}$ and the source distribution $I_{S}(x,y)$ on $z=z_{0}$ can then be calculated numerically, e.g. from the law of refraction and the Jacobian equation (2). $I_{S}(x,y)$ is shown in Fig. 9(b) and serves as the input for the calculation of $\mathbf{u}^{\infty}(x,y)$ ( $417sec$ ).

For the calculation of the intial surface with Eq. (15), we have to fix one point of the freeform surface relative to the plane $z=z_{0}$ , which we choose to be $z(0,0)=z_{0}+1$ and the side length of the target area is chosen to be $2$ . The calculation of $z_{S}^{\infty}(x,y)$ took $145sec$ and the root finding of Eq. (10) about $215sec$ . Since the input vector field is noncollimated, the surface $z(x,y)$ is again calculated on a scattered grid and has to be interpolated for the ray-tracing simulation. The results of the ray tracing are shown in Fig. 10 and Table 3. Also for this example, the results show that the initial map is appropriate for the root finding of Eq. (10) and that by the presented design algorithm, single freeform surfaces for complex target irradiances and nonplanar or nonspherical input wavefronts can be calculated.

The quantitative differences between the ray-tracing results for the examples are thereby mainly due to interpolation and extrapolation effects. It is also a priori not clear, if for the particular predefined system parameters of each example an exact solution to Eq. (10) exists. This might also lead to a degradation of the qualitiy of the ray-tracing results, since an accurate root-finding of the respective nonlinear equation systems might not be possible. In contrast to the first two examples, additional errors are introduced in the third design example by the numerical calculation of $I_{S}(x,y)$ and the input direction vector field.

VII Conclusion

A single freeform surface design algorithm for prescribed input wavefronts and complex irradiance distributions was introduced. This was done by deriving the PDE system (10) and Monge-Ampère equation (11), which governs the single freeform design process for ideal sources in the geometrical optics limit, from the energy conservation and ray tracing equations in cartesian coordinates. After introducing a possible way to calculate an appropriate initial iterate in section IV, the PDE system is then solved for the unknown functions $\mathbf{u}(x,y)$ and $z(x,y)$ by a finite difference discretization and a root finding of the resulting nonlinear equation system. The design examples in section VI show that the presented design algorithm extends the results of Ref. Boe16_1 for collimated beam shaping to the nonparaxial regime and that the algorithm is applicable to point sources and more general input wavefronts. This might be helpful in applications in which one has to deal with the shaping of e.g. astigmatic wavefronts produced by laser diodes or simplify manufacturing processes. Also the relaxation of physical situations, e.g. for point sources with large opening angles, might be reached by a predefined entrance surface of a freeform lens.

Despite the presented results, there are still interesting problems left for future research and improvements, respectively. For example, a mathematical rigorous justification of the initial iterate construction method presented in section IV is of interest, since there is no guarantee that the proposed method will always give a suitable intial iterate, despite the arguments given in section IV. Also the numerical extension of the presented algorithm to more complicated boundaries of the irradiance distributions would help to increase the number of possible applications or energy efficiency, respectively. This requires the numerical implementation of the transport boundary conditions, which are realized for irradiance distributions with square boundary shapes by Eq. (19), for more complex boundary shapes and will be investigated in the future. As mentioned in section III, a numerical investigation of the Monge-Ampère equation in Eq. (11) is of interest. Despite the listed additional complications, the reduced number of design variables can benefit the speed of the root-finding process. Another interesting future goal is the generalization of the presented design algorithm to wavefront and irradiance control with double freeform surfaces by combining it with the results presented in Ref. Boe17_1 and the investigation of the manufacturability of the calculated freeform surfaces.

Funding

We acknowledge support from the German Federal Ministry of Education and Research (BMBF) under contract number 031PT609X within the project KoSimO.

Acknowledgments

The authors want to thank Norman Worku for providing and helping them with the MatLightTracer and Johannes Stock for comments on the manuscript.

Bibliography21

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1(1) D. Michaelis, D. Schreiber, and A. Bräuer, “Cartesian oval representation of freeform optics in illumination systems,” Opt. Lett. 36 (6), 918–920 (2011).
2(2) V.I. Oliker, “Mathematical aspects of design of beam shaping surfaces in geometrical optics,” in Trends in Nonlinear Analysis , M. Kirkilionis, S. Kromker, R. Rannacher, and F. Tomi, eds. (Springer-Verlag, 2003), pp. 193-222.
3(3) F.R. Fournier, W.J. Cassarly, and J.P. Rolland, “Fast freeform reflector generation using source-target maps,” Opt. Express 18 (5), 5295–5304 (2010).
4(4) V.I. Oliker, “Controlling light with freeform multifocal lens designed with supporting quadric method(SQM),” Opt. Express 25 (4), A 58-A 72 (2017)
5(5) H. Ries and J. Muschaweck, “Tailored freeform optical surfaces,” J. Opt. Soc. Am. A 19 (3), 590–595 (2002).
6(6) R. Wu, L. Xu, P. Liu, Y. Zhang, Z. Zheng, H. Li, and X. Liu, “Freeform illuminaton design: a nonlinear boundary problem for the elliptic Monge-Ampère equation,” Opt. Lett. 38 (2), 229–231 (2013).
7(7) K. Brix, Y. Hafizogullari, and A. Platen, “Designing illumination lenses and mirrors by the numerical solution of Monge–Ampère equations,” J. Opt. Soc. Am. A 32 (11), 2227-2236 (2015)
8(8) R. Wu, P. Liu, Y. Zhang, Z. Zheng, H. Li, and X. Liu, “A mathematical model of the single freeform surface design for collimated beam shaping,” Opt. Express 21 (18), 20974–20989 (2013).