Design of a freeform two-reflector system to collimate and shape a point source distribution

In this paper we propose a method to compute a freeform reflector system for collimating and shaping a beam from a point source. We construct these reflectors such that the radiant intensity of the source is converted into a desired target. An important generalization in our approach compared to previous research is that the output beam can be in an arbitrary direction. The design problem is approached by using a generalized Monge-Amp\`ere equation. This equation is solved using a least-squares algorithm for non-quadratic cost functions. This algorithm calculates the optical map, from which we can then compute the surfaces. We test our algorithm on two cases. First we consider a uniform source and target distribution. Next, we use the model of a laser diode light source and a ring-shaped target distribution.


Introduction
Beam shaping is an important research topic within illumination optics. Especially the shaping of a parallel beam into another parallel beam with a different distribution is well researched [1][2][3]. This is often linked to the shaping of laser beams. A common source for a laser is a laser diode, which can be modelled as a point source [4,5]. The diverging beam from such a diode is often collimated first with a lens, before other manipulations such as beam shaping are applied [5]. An optical system that directly shapes and collimates the output from a point source is able to skip the first collimation step. This can reduce the total number of necessary optical surfaces and increase the efficacy by avoiding Fresnel reflections. An example of a useful light distribution in a collimated beam that we will discuss is the ring shape, meaning that the projection of this beam on a plane perpendicular to it gives a ring-shaped illumination pattern. A possible benefit of such a ring-shaped target is that a subsequent focussing can be done more accurately [6]. A ring-shaped illumination pattern also has a use in welding [7].
In most research regarding collimated beams, it is assumed that the outgoing beam is in the same direction as the incoming beam. For the case of a point source, however, there is not one single direction of light emission, but often there exists a symmetry axis that plays the role of beam direction. In this paper we will drop that assumption and allow the outgoing beam to be in any arbitrary direction. This gives us the possibility to create so-called folded optics. As a result we can design more compact optical systems. As mentioned, we will look for an optical system to collimate a beam from a point source. For this we need two optical surfaces, one to shape the light to the desired target distribution and one to collimate the beam. We choose to work with two reflector surfaces as freeform optical surfaces.
Optical design for illumination can roughly be divided into two categories: forward and inverse methods. The former deals with the calculation of the output of an optical system. The result of such a forward method can then be used to iteratively refine the design, which is a slow process [8]. In this paper we will focus on inverse methods. With such methods, the goal is to compute the shape of the optical system given source and target illumination patterns.
Although the example of a laser gives an idea of the possibilities, we will not restrict ourselves to this. In our derivations we will use the approximation of geometrical optics. This way we can view the calculation of the surfaces as an optimal transport problem, i.e., to find an optical map that 'transports' the light from the source to the target. This gives us a Monge-Ampère type equation with transport boundary condition [2].
Several methods have been developed for solving problems similar to the point source and parallel outgoing beam. A more thorough overview of inverse methods can be found in [9]. Some of them use numerical methods such as finite differences and Newton's method to directly solve the Monge-Ampère equation [10,11]. Oliker et al. proposed the supporting quadric method [1,12]. Alternatively, the optimal mass transport problem is reduced to a linear assignment problem by Doskolovich et al. [13]. Another approach by Feng et al. uses ray mapping to calculate the shapes and positions of the surfaces [14]. A ray mapping method has been used to construct optical systems with two freeform surfaces for arbitrary input and output wavefronts [15].
To the best of our knowledge, none of the above approaches have been used to specifically design an optical system for collimating and shaping a diverging beam from a point source. In this paper we will modify a least-squares method to solve the Monge-Ampère type equation. Versions of this algorithm have been used before for multiple optical design challenges, including parallel to far-field [16], parallel to parallel [17] and point source to far-field [18].
In Section 2 we present the mathematical model linking the shapes of the surfaces to the source and target distributions. The equations in this model are used in the algorithm mentioned before. We will give a short summary of this algorithm in Section 3. For some parts of the algorithm we refer to other papers and we briefly discuss the most important parts. In Section 4 we test our algorithm on two test cases: First a uniform point source and a uniform target, second a laser diode source and a ring-shaped target. The conclusion of our findings is given in Section 5.

Formulation of the mathematical model
In this section we formulate the mathematical model of a reflective optical system creating a collimated beam from a point source. Given a source light distribution we want to design an optical system of two reflectors. The first reflector will be used to shape the intensity profile and the second one will collimate the beam. The collimated output beam should give a light distribution on a target plane at a distance l from the source, where the distribution is a given function of the position coordinates in the plane. A two-dimensional illustration of the system is given in Fig. 1.

Derivation of the cost function
We choose the point source to be located at the origin of a coordinate system given by the standard basisê 1 ,ê 2 ,ê 3 . The orientation for this coordinate system is arbitrary. Often, thê e 3 -axis is chosen to coincide with the symmetry axis of the bundle of light emitted from the source. However, such a symmetry does not necessarily exist and will not be assumed for our derivations. The only restriction we have on the coordinate system is that there can be no light emitted in the negativeê 3 -direction. The reason for this will become apparent later. This source emits rays with unit direction vectorsŝ = (s 1 , s 2 , s 3 ) T . Two reflectors, R 1 and R 2 , are used to collimate and shape the light emitted by the source into an output beam with the desired intensity distribution. These output rays should propagate as a parallel beam that is parallel to a given direction vectorâ 3 perpendicular to the target plane. The direction ofâ 3 with respect toê 1 , . . . ,ê 3 is given by a rotation ofê 3 with a polar angle ϕ around theê 2 -axis and subsequently with an azimuthal angle θ around theê 3 -axis. The matrix associated with this composite rotation is given by and we defineâ i = Aê i , i = 1, 2, 3. Theseâ i -vectors are the rotations of an orthogonal basis, so they are again an orthogonal basis andâ i is equal to the ith column of A. The target plane is perpendicular toâ 3 at a distance l from the point source alongâ 3 . Any point on the target plane is then given by y 1â1 + y 2â2 + lâ 3 . The position in the target plane is denoted by y = (y 1 , y 2 ) T . The first reflector, R 1 , is defined by the radial distance from the source, u = u(ŝ). The second reflector, R 2 , is given by the perpendicular distance from the target plane, w = w(y). The reflectors are then mathematically described by where r 2 (y) = y 1 , y 2 , l − w(y) T . With P and Q we denote the points where a ray hits the first and second reflector, respectively, see Fig. 1. The distance between those two points is denoted by d. In the framework of optimal mass transport we want to derive a relation betweenŝ and y of the form where u 1 and u 2 are related to the shape and location of the surfaces and c is the so-called cost function [2]. For any ray, we introduce the following notation. The 2-vectors q s = 0 and q t = y denote the position of the ray at the source and target, respectively. These are projections of 3-dimensional position vectors onto the plane z = 0 (source plane) and the target plane, respectively. Similarly, p s and p t = 0 are the projections of the direction vectors on the source and target planes. In terms of Hamiltonian characteristics, the optical path length (OPL), denoted by L, is equal to the point characteristic V [19]. We can write Because the variablesŝ and y denote a direction at the source and a position on the target, we work with the second mixed characteristic W * = W * (p s , q t ), given by However, since q s = 0 we can dismiss the second term and the mixed characteristic is equal to the optical path length. The following relations can be derived for the mixed characteristic [19] This proves that W * is independent of the direction from the source and the position on the target plane. As a result, the optical path length L = V = W * is a constant. We will eliminate d from Eq. (4) by using the fact that d is the distance between P and Q. We denote by r P = r 1 (ŝ) and r Q = r 2 (y) the position vectors of P and Q, respectively. Furthermore, we use r Q = r 2 (y), so we can write Because A is a rotation matrix, we have |Ar Q | = |r Q |. We also have r P = uŝ and A = (â 1 ,â 2 ,â 3 ). This can be used to write where A 2 = (â 1 ,â 2 ) ∈ R 3×2 and we omit the dependence of u and w onŝ and y for now. We combine Eq. (4) and Eq. (8) to obtain We want to separate u and w to get an equation of the form of Eq. (3). For that, we first divide Eq. (9) by u. It is reasonable to assume u > 0, since otherwise the reflector coincides with the source at some point. Introducingũ = 1/u, Eq. (9) reads The parameter β = L − l is introduced, which is called the reduced optical path length. This is used to rewrite the equation above as To obtain an equation of the desired form, we will factorize the equation above and then apply logarithms to both sides of the equation. The factorization leads to As mentioned, we would like to take the logarithm of both sides of the equation to get to the form of Eq. (3). However, we need to make sure that both sides are positive. We use L = u + d + w and we can write wheret is the (unit) direction vector after the first reflection. Combining both relations, we Sincet is a unit direction vector, we know that β ≥ (1 −ŝ ·â 3 )u ≥ 0, with the first equality only ift ·â 3 = 1. In that caset is parallel toâ 3 , so there would be no second reflection. We can disregard this case and therefore write The vector y is the displacement along the plane spanned byâ 1 andâ 2 . We have y 1â1 +y 2â2 = A 2 y. This displacement is determined completely by the first two ray segments; from the source to the first reflector, and from the first to the second reflector. Using this, we write Substituting this into the expression for κ 2 and using Eq. (4) and Eq. (13) gives us We now use again that (17) Note that we can combine theŝ ·â i terms into a vector written as This vector can be used to eliminate the u 2 term, since as A is a rotation matrix and therefore AA T = I. Similarly, swappingt forŝ gives us 3 i=1 (t ·â i ) 2 = 1. We use this to rewrite with equality if and only ift =ŝ. This will not occur, since there would be no reflection at the first surface in that case. We deduced that κ 1 > 0 and κ 2 < 0. Consequently, T must be negative as well. Both sides of Eq. (12a) are therefore multiplied by −1 to obtain Before we apply the logarithm to this equation, we scale all the lengths by a factor β. Note thatŝ is already dimensionless. We introduce the variable z such that y = βz. We substitute this into the function −T (ŝ, y) =T (ŝ, z) and obtain Next, we scale all the other lengths, viz. w(y) = βŵ(z), L = βL and l = βl. Furthermore,ũ has the dimension of length inverse, so we scale this by 1/β toũ(ŝ) = (1/β)û(x). Substituting this in the expressions for κ 1 and κ 2 gives So we can now define the new functionsκ 1 andκ 2 bŷ and For the algorithm that will be introduced in Section 3, we want to change the source coordinates into two independent variables instead of three variables on the unit sphere. For that we choose the stereographic projection from the south pole (0, 0, −1) onto the equator plane z = 0, written as x. These coordinates are given by Furthermore, we introduce When we take the logarithm of both sides of Eq. (25), taking into account the functions we just introduced, we get the desired form of Eq. (3), given by where the function c is called the cost function and is defined as

Energy conservation
We have deduced an equation that implicitly links x and z. We assume that there exists an explicit relation given by z = m(x). A constraint for this function is that the energy of the source should be conserved through the optical system. By S ⊂ S 2 , with S 2 the unit sphere, we denote the set of direction vectorsŝ from the source. For the stereographic projection x of this set we use X ⊂ R 2 . The energy density of the source is given by f = f (ŝ). We write T ⊂ R 2 for the set of target positions y and Z ⊂ R 2 for the set of scaled target variables z. The desired illuminance on the target plane is then given by g = g(y) for y ∈ T . The energy of any subset of the source domain should be conserved through the optical system. For any subsetÃ ⊂ S there is a corresponding set A of stereographic coordinates. Energy conservation is given by the equation where y = βz. The first integral is over a surface element on the unit sphere, while the second one is over an area element in R 2 . When A = S we have the special case of global energy conservation. This means that the total energy in the source and target should be equal. We assume that f and g are constructed such that this is true. First, we want to transform the left-hand side integral to an integral over x instead ofŝ. From integration by substitution we have Now, for the right-hand side of Eq. (30) we also want to transform the integration variable to x. For that we first have to change from y to z, since we know z = m(x). We had defined y = βz, so the Jacobi matrix of this transformation is given by β 2 . Substituting this and Eq.
We use Dm to denote the Jacobi matrix of m. This equation should hold for any A corresponding to anÃ ⊂ S. Therefore the integrands should be equal almost everywhere. Then for any x(ŝ) withŝ ∈ S we have the generalized Monge-Ampère equation The boundary condition of the problem is a transport boundary condition given by It states that the boundary of the source domain should be mapped to the boundary of the target domain [20]. For the remainder of this article we will assume that det(Dm) > 0, so we can ignore the absolute value in Eq. (33a).

Polar coordinates
In general, the rays will be emitted by the source in a conical bundle symmetric around the z-axis. In that case, the source domain in stereographic coordinates, X , will be a circle. So, it makes sense to switch to the polar coordinate system. The polar stereographic coordinates are written as ω = (ρ, ζ) and the transformation x = x(ω) is given by We define Ω to be the source domain in polar stereographic coordinates, so Ω = ω(X ). Furthermore we define u * 1 (ω) = u 1 (x(ω)) and c * (ω, z) = c(x(ω), z), so equation (28) changes to u * 1 (ω) + u 2 (z) = c * (ω, z).
We also need to transform the energy conservation equation to polar coordinates. With integration by substitution, we change Eq. (32) to where dω = dρdζ. We introduce the notation m * (ω) = m (x(ω)), implying z = m * (ω). Note that Eq. (36) should hold for A corresponding to any arbitrary part of the source domain, so for any ω ∈ Ω we have The matrix Dm * is the Jacobian of m expressed in polar coordinates and is given by This follows from a coordinate transformation of the Jacobian with Cartesian coordinates in Section 2.2. The boundary condition (33b) changes to ∂Z = m * (ω (∂X )) .
Note that we have Ω = ω(X ), but ∂Ω = ω (∂X ). For example, the relation ρ = 0 would give a point in the interior of X , but a line on the boundary of Ω. In the remainder of this article we will omit the asterisks, because we will only work with functions of ω.

Mapping
We want to find an (implicit) expression for the optical mapping z = m(ω). This can be done using Eq. (35). There are many solutions to solve this equation for m, so we make a special choice [2]. It is possible to find a c-convex pair of functions u 1 , u 2 such that Conversely it is possible to find a c-concave pair defined by In either case, the solution for u 2 has an argumentω that is a stationary point. This leads to the requirement that where the gradient with respect to ω is given by To ensure that the solution is a maximum we need the Hessian w.r.t. ω of c(ω, z) − u 1 (ω) to be symmetric negative definite (SND). Similarly, to ensure a minimum we require the Hessian to be symmetric positive definite (SPD). The Hessian matrix of any function v in polar coordinates is given by [21] H In the Hessian of the function c(ω, z) − u 1 (ω) the first derivative terms cancel because of Eq.
(40), so this Hessian is given by Note that D ωω v is not equal to the Hessian matrix of v. We assume that a mapping z = m(ω) exists, and substitute this into Eq. (40). We then take the derivative with respect to ω and apply the chain rule to obtain where (45b) Note that the first term in Eq. (45a) means differentiating c(ω, z) w.r.t. ω twice and then substituting z = m(ω). We can rewrite Eq. (45a) to Note that −P is the Hessian matrix in Eq. (43). Therefore, for a c-convex or c-concave pair of functions u 1 , u 2 , we have the condition that P should be SPD or SND, respectively. We assumed that a mapping m is defined (implicitly) by Eq. (40). By the implicit function theorem such a mapping is guaranteed to exist if the Jacobian matrix of the left-hand side of the equation with respect to z is invertible [22,Sec. 12.8].
To summarize, we need to find a mapping m satisfying Eq. (37c) and Eq. (46) for a matrix P that is SPD or SND, with det(P ) = F det(C). The exact mapping is determined by the parameters l and β (or L) and the direction of the outgoing beam, given by ϕ and θ.

Numerical method
We explain the least-squares algorithm which we use to solve the problem derived in the previous section. This algorithm has been explained thoroughly for a source with Cartesian [2,16] and polar [9] coordinates. In this section we give a brief overview. We first compute the mapping, followed by a calculation of the surfaces. We restrict ourselves to the c-convex solution of Eq. (35). As shown in the previous section, we need to solve where C = D ωz c and P satisfies det(P ) = F det(C), see Eq. (37a). To get the c-convex solution, P needs to be SPD. We enforce the equality in Eq (47) by minimizing the functional J I defined by The norm . F is the Frobenius norm. To enforce the boundary condition we minimize the difference between the (given) boundary of T and the mapping of the boundary of the source: where b ∈ ∂Z. We combine these two functionals into a weighted average with parameter α ∈ [0, 1] given by All of these functionals are defined on the following spaces We cover the domain Ω by a grid with gridpoints ω ij . The algorithm to find m is initialized by a guess m 0 for the mapping. With this mapping we compute the matrix C 0 . Then, we iteratively perform the next steps either for a fixed number of iterations or until a stopping criterion is met, The minimization procedures for b and P do not contain derivatives of their respective variables, so these can be minimized pointwise. A method for solving step (52a) is given by Romijn et al. [18]. We will explain step (52b) in a bit more detail. The matrix P needs to be SPD. The symmetry of P is enforced by defining P = p 11 p 12 p 12 p 22 .
We approximate Dm using central differences and define Q = CDm. Instead of minimizing Q − P F we solve an equivalent problem with the same minimizers [16]. We introduce the symmetric matrix Q S = 1 2 (Q + Q T ), with off-diagonal entries q S = 1 2 (q 12 + q 21 ). The optimization problem that needs to be solved is then minimize H S (p 11 , p 22 , p 12 ) = 1 2 Q S − P 2 F , subject to det(P ) = F det(C).
It turns out that we can always select at least one solution of this problem that satisfies the constraint that P is SPD [2]. To solve problem (54) we use the Lagrange multiplier method. We introduce the Lagrangian function Λ(p 11 , p 22 , p 12 , λ) = H S (p 11 , p 22 , p 12 ) + λ p 11 To find stationary points we take the derivatives w.r.t. each variable and set them equal to zero. Elementary calculation lead us then to the system of equations Solutions of this system can be calculated analytically and explicitly [2]. In the case that we find multiple solutions, we have to choose the one that gives the lowest value for H S .

Computing the mapping
In the functional J there are derivatives of m, so we can no longer optimize pointwise. To be able to minimize this functional, we apply calculus of variations. The first variation of m in the direction of an arbitrary function η = (η 1 , η 2 ) T ∈ M is given by where : denotes the Frobenius inner product associated with the Frobenius norm . F . Like before, Dη denotes the Jacobian of η w.r.t. ω. We can rewrite the inner product in the first integral to obtain (CDm − P ) : where w 1 and w 2 are defined by We use this and Gauss's divergence theorem to rewrite the integral over Ω as wheren is the outward unit normal of X . We substitute this into Eq. (57) and set the first variation equal to zero to obtain At a minimum of J, the first variation should be equal to zero for any vector η. We can split this in two cases. First we choose η 2 = 0 and set the first variation equal to zero for any η 1 ∈ C 2 (Ω). Similarly, we have the case where η 1 = 0. With the use of the fundamental lemma of calculus of variations [23] we get the boundary value problem div(C T CDm) = div(C T P ), for ω ∈ Ω, (62a) where div is defined in the following way. Let B = (b ij ) ∈ R 2×2 , then div(B) = 1 ρ (63) This boundary value problem is then solved using the finite volume method [9, App. A].

Computation of the reflector surfaces
The algorithm (52) computes a mapping m. From this mapping we can compute the shape of the surfaces, given by u and w. To find u we first compute u 1 from equation (40). We introduce the functional I to quantify how close a function φ is to the exact solution of that equation. Our solution u 1 is then given by where To solve this optimization problem we use calculus of variations. The first variation of (65) is analogous to Eq. (57). Similar to the calculation of the mapping, we use Gauss's divergence theorem and the fundamental lemma of calculus of variations to get the boundary value problem ∆u 1 = div(∇ ω c), ω ∈ Ω (67a) Here, ∆ denotes the Laplace operator, which in polar coordinates is given by Let u 1 be a solution to this boundary value problem. To calculate the shapes of the reflector surfaces, u and w, we combine Eq. (24) and Eq. (27) to obtain where u 2 = c − u 1 , using Eq. (28). The solution to the boundary value problem (67) is unique up to an additive constant [9]. This means that choosing the value of u 1 in one point gives us a unique solution. We use this degree of freedom to choose the position of one of the reflectors along the central ray given byŝ =ŝ 0 := (0, 0, 1) T . For example, we can choose u(ŝ 0 ) = u 0 for some u 0 > 0. With Eq. (69), this then gives where ω 0 = ω(ŝ 0 ) and ϕ is the polar angle of the output rays. This expression, together with the boundary value problem (67), gives a unique solution for u 1 and therefore also for u and w. In this case we cannot freely choose the position of the second reflector along the central ray. Alternatively, we could give w m(ω 0 ) = w 0 as input instead of u 0 . This would also fix u 1 , u and w.

Results
We will apply the method from the previous section to examples with two distinct combinations of source and target distributions. We create several optical systems with different layouts to show the possibilities of the algorithm.

Uniform to uniform
The example we will look at first consists of a uniform conical source and a uniform circular target distribution. We will test several layouts of optical systems, so different angles for the outgoing beam. Because our source is rotationally symmetric, our choice of θ does not matter and we can choose θ = 0. We will discuss optical systems with ϕ = 0, ϕ = π/2 and ϕ = π. The first and last one have outgoing beams parallel to theê 3 -axis in the positive and negative direction, respectively. For ϕ = π/2 the output beam is parallel to theê 1 -axis. We shift the target domain along the target plane to avoid mirrors obstructing rays. Otherwise, for example, the first reflector could be in the way of the rays from the second reflector to the target. For the point source, we use a uniform distribution in the direction vectorŝ rather than the stereographic coordinates x. This means that the intensity distribution given by f = f (ŝ) is constant. The source domain is given by ρ ≤ 0.1 and we choose the value of f such that the total flux is 1. The target domain is a circle with radius 2. The target distribution g is constant over this domain, with a flux of 1.
There are some parameter choices that will affect the optimization procedure or the resulting optical system. For example, the functional in equation (50) contains the parameter α. The smaller α, the more important the boundary is in the optimization, relative to the interior. The choice of α will have an impact on the convergence speed of the algorithm [2]. Other parameters have an influence on the layout of the optical system. These are the distance to the target plane, l, the reduced optical path length, β (or L), and either u 0 or w 0 .
For our first test case we use α = 0.1. From numerical experiments we found that this value works well. We discretize the polar source domain with a 200 × 200 grid. The results of our test case were obtained by running the least-squares algorithm for 200 iterations. First, we did this for a case where ϕ = 0. We use a target plane at a distance l = 20 and we choose to shift the target domain by a distance −10 along theâ 1 -direction. Furthermore, we put u 0 = 10 and β = 15. The resulting optical system is visualized using a ray trace procedure. The result can be seen in Fig. 2. In the figure we show 100 random rays that are obtained from this ray trace. On the target plane we show the illumination pattern in a bounding box of the target domain. We divide this box in 100 by 100 bins and use a quasi-Monte Carlo method tracing ten million rays to get an illumination pattern. A better view of the illumination pattern is given in Fig. 3. As we can see, the illuminated part of the bounding box forms a circle. The flux per bin, and thus the intensity, on this circle is nicely uniform. The bins within the circle have some variation in flux that is caused by the Monte Carlo method.
As mentioned before, an important benefit of our algorithm is the ability to construct a wide variety of optical systems. Some examples of possible layouts are shown in Fig. 4. The first two figures both have ϕ = π/2, θ = 0 and β = 22. The target plane is located at a distance l = 5 and the center of the target domain is shifted by −5 alongâ 1 . With ϕ = π/2 and θ = 0 we haveâ 1 = −ê 3 . The difference between the two optical systems is that in 4a we have u 0 = 10 and in Fig. 4b we have u 0 = 1. Note that these two have the exact same mapping, since the parameter u 0 only plays a role in the calculation of the surfaces (see Sec. 3). Therefore, we only need to calculate the mapping once, and we are able to construct both of these optical systems from this mapping by varying u 0 . For Fig. 4c we change the angle of the outgoing beam to ϕ = π and we use w 0 = 6. We choose l = 0, so the target plane is equal to the source plane. We shift the target domain by 10 along theê 1 -axis and we again use β = 22. In each of the previously mentioned optical systems, there is a plane of symmetry. However, this is not necessary for our algorithm. In Fig. 4d, we shift the target domain by −5 alongâ 1 and 5 alongâ 2 . This breaks the symmetry of the optical system. Furthermore, we set β = 15 and w 0 = 5.

Laser diode to ring-shaped target
The second test case for our algorithm is the case of a laser diode [5] to a ring-shaped target pattern. In this case we choose ϕ = π/2 and θ = 0, so that the outgoing beam is parallel to thê e 1 -axis. Note that the angle θ only matters when the source distribution is not rotationally symmetric. The intensity distribution of a laser diode can be modeled by an elliptical Gaussian on a plane perpendicular to theê 3 -axis [4]. The emitted light has a 1/e 2 intensity angle of θ x in theê 1 ,ê 3 -plane, and θ y in theê 2 ,ê 3 -plane. From this we can derive a density in terms of the x-variables, see App. A. The following density function is obtained: with σ x = tan(θ x /2) and σ y = tan(θ y /2). The scalar B x is a scaling parameter such that the flux of f over X is 1. The parameter δ is used to ensure a minimal value for the source density. The source domain X is given by where Θ = max(θ x , θ y )/2. We want to have a target which consists of a ring on a target plane with an outer radius r o and an inner radius r i . For this we will construct a density as a function of the unscaled  target variables y. We want this density to be uniform on a ring Ω r . With our algorithm we need a simply connected domain. So, instead of a ring, the domain we use will be a disk Ω T = Ω r ∪ Ω i , where Ω i is the circular domain enclosed by Ω r . The ring will have a higher density than the inner disk. Our model assumes that the density is smooth, so we will need to approximate this discontinuous change in density between the ring and the inner circle. From the derivation in App. B we obtain the target intensity distribution We choose c so that the total flux over the target domain is equal to 1. The parameters k and ε determine how close this density is to the ideal density, uniform on Ω r and zero on Ω i . They also influence the convergence of the algorithm. Generally a smaller k or larger ε will give better convergence, but a less pronounced difference in density between the inner region and the outer ring.
For the source density, we choose to work with 1/e 2 intensity angles θ x = 45 • and θ y = 13 • . These are typical values for a laser diode [5]. To avoid a too large difference in intensity in the source, we choose δ = 10 −3 in the density function from Eq. (71). The target domain consists of an outer radius r o = 2 and an inner radius r i = 1. This is shifted by −5 along theâ 1 -axis in the target plane at a distance l = 5 from the point source. We choose ε such that the density in the inner circle is at least 10% of the density in the ring. Furthermore, we used k = 100 in the target density. We set β = 30 and u 0 = 20. Because the intensity on the boundary of the source is much lower than at the center, it is more difficult for the algorithm to find a good mapping for the boundary compared to our previous test case. To counter this problem we decrease α. We now choose α = 0.01 to put more emphasis on J B relative to J I . Compared to the previous cases, we also increase the number of grid points and apply a 400 × 400 grid to the source domain. Tracing one million rays through the system that results from our algorithm gives us the illumination pattern and optical system in Fig.  5 and Fig. 6.

Conclusion
In this paper we introduced a method for computing the shapes of two reflectors to collimate a beam from a point source, for given source and target light distributions. A specific point of interest is the fact that the outgoing beam can be in any arbitrary direction. First we have derived a relation between the surface shapes and locations, and the optical mapping. Then we derived a Monge-Ampère type equation for this problem. We proposed an algorithm based on the least-squares method to solve this equation.
We tested our algorithm on two test cases. A uniform source and target, and a model of a laser diode to a ring-shaped target. The former was used to test the algorithm and the different layouts that it could attain. The latter consisted of more complicated source and target distributions instead. The second test case showed us that we might need to put some restrictions on our model to improve the convergence of the algorithm. For example, the variation of the intensity in the source cannot be too large. This was solved by imposing a minimal value of the source distribution.
For further research there are several points of interest. In this paper we have skipped over some practical constraints. There is nothing in the algorithm yet preventing the rays from crossing a reflector. For example, the first reflector might (partially) be in the way of the outgoing beam. We managed to work around this by varying parameters to get a feasible solution. However, for the future it might be interesting to research if it could be possible to incorporate physical constraints like this into our algorithm. Furthermore we want to extend our algorithm to include more physical phenomena. Examples of this are Fresnel reflection (for lenses) and scattering.

B RING-TARGET DENSITY
We write the right-hand side as an integral over A by using integration by substitution. Because this holds for any subset A, this gives a source density function f , with f (x) = det Dp(x) Ĩ p(x) This density seems to be dependent on z p , while the intensity of the point source should of course not depend on the plane of projection. We will show that we can in fact write f as a function independent of z p . In general, the increase of w 1 and w 2 for increasing z p is given by the full angles at the source point, θ x and θ y . So these radii of a laser diode increase linearly in z p . We then have w 1 = z p tan(θ x /2), w 2 = z p tan(θ y /2), so that with σ x = tan(θ x /2), σ y = tan(θ y /2) and B x = 4Bz 2 p . With Eq. (75) and integration by substitution, B x can be written as This shows that B x is independent of z p and indeed also our source density f is independent of z p . The only part left for our source density is to define the boundary of the domain in stereographic coordinates. We mentioned before that the boundary of the domain in pcoordinates on the plane z = z p is given by |p| = max w 1 , w 2 , which is written as |p| = max z p tan(θ x /2), z p tan(θ y /2) = z p tan (Θ/2), with Θ = max θ x , θ y /2. From Eq. (76) we obtain |x| = z 2 p + |p| 2 − z p |p| .
Again, this result is independent of the choice of z p . We have now modeled the laser diode by the intensity distribution in Eq. (80), defined on the domain X = x ∈ R 2 | |x| ≤ tan (Θ/2) .

B Ring-target density
We define our circular domain Ω T by an angle ξ ∈ [0, 2π) and a radius r ∈ [0, r o ]. The outer ring is denoted by Ω r and given by r ∈ [r i , r o ]. The inner disk Ω i is then given by r ∈ [0, r i ] such that we have Ω T = Ω r ∪ Ω i . We denote the area of the ring by A and we have