Numerically solving generated Jacobian equations in freeform optical design

We present an efficient numerical algorithm that can be used to solve the generalized Monge-Ampère equations for a single freeform reflector and lens surface. These equations are instances of so-called ‘generated Jacobian equations’ which are characterized by associated generating functions. The algorithm has a wide applicability to any optical system that can be described by a smooth generating function.


Introduction
We consider two optical systems: (System 1) a lens surface for a point source and far-field target intensity, and (System 2) a reflector surface for a parallel incoming beam and near-field target illuminance. For System 1 the source emits light radially outward and we consider a source emittance in spherical coordinates f 1 (φ, θ) in [lm/sr]. The target intensity distribution g 1 (ψ, χ) in [lm/sr] in the far field is expressed with respect to a different set of spherical coordinates (ψ, χ) with origin the lens surface approximated as a point in space (far-field approximation). The first surface of the lens is spherical and the second surface is freeform. For System 2 we are given an emittance of the source domain f 2 (x) in [lm/m 2 ] and a target illuminance in the near field given by g 2 (y) in [lm/m 2 ], where x and y are the local Cartesian coordinates of the source and target planes, respectively. In figure 1 a schematic representation is given for System 2. Our goal is to compute the optical surfaces that define the mappings m 1 and m 2 transforming f 1 into g 1 and f 2 into g 2 , respectively.

Mathematical model
In previous work on System 1 [1], we used optimal transport theory to derive a relation of the form u 1 (x) + u 2 (y) = c(x, y), where u 1 specifies the location of the freeform optical surface and c(x, y) is a logarithmic cost function in optimal transport theory with (x, y) denoting the stereographic coordinates obtained from performing coordinate transformations on the source and target domains. Subsequently, we could derive an optical map m and derive the generalized Monge-Ampère equation as 2 , * e-mail: l.b.romijn@tue.nl * * e-mail: m.j.h.anthonissen@tue.nl * * * e-mail: j.h.m.tenthijeboonkkamp@tue.nl * * * * e-mail: wilbert.ijzerman@signify.com where Dm is the Jacobi matrix of m, P is a symmetric positive-definite (SPD) matrix and C = D xy c = (c x i y j ) is the matrix of mixed second-order partial derivatives with respect to the stereographic coordinates of the source and target. Figure 1: Schematic representation of a parallel beam of light rays, reflected by the reflector surface onto a target in the near field.
However, for some optical systems such as System 2 a relation of the form u 1 (x) + u 2 (y) = c(x, y) does not exist in any coordinate system. In this paper, we generalize the cost-function formulation in optimal transport theory to a generating-function formulation by defining u(x) = G(x, y, z(y)) and z(y) = H(x, y, u(x)), where u is the unknown location of the optical surface, G and H are smooth functions and G(x, y, ·) and H(x, y, ·) are each other's inverses. For both systems, we can reformulate the generalized Monge-Ampère equation as a generated Jaco- bian equation [2] det where P is a symmetric positive-definite (SPD) matrix and we redefine C as C = D xy H(x, m(x), u(x)) = (H x i y j ). Note the added dependency of C on the surface u. For System 1, we consider the stereographic coordinates (x, y) and have F(x) = f 1 (x)/(1 + |x| 2 ) 2 and G(m(x)) = g 1 (m(x))/(1 + |m(x)| 2 ) 2 . For System 2, we take Cartesian coordinates x = x and y = y, and derive F(x) = f 2 (x) and G(m(x)) = g 2 (m(x)).

Numerical method
Using a generalized least-squares approach [4] we iteratively update the mapping m. We write the generalized Monge-Ampère equation as the matrix equation CDm(x) = P(x), where P(x) satisfies det(P(x)) = F(x)/G(m(x)), cf. [5]. An iterative procedure is used to find the mapping, which involves an efficient procedure to find the numerical solution of a constrained minimization problem for each grid point to compute P and the solution of a linear elliptic boundary value problem to compute m.
In order to update the matrix C during the iterative procedure, the location of the optical surface is calculated from the mapping also in a least-squares sense. The numerical method can be used to compute reflector and lens surfaces for various examples. For instance, we can compute the reflector surface that transforms light of a parallel beam into a near-field picture of Porto, as shown in figure 2. For System 1 we will compare the results of the cost-function approach previously presented in [1] to the new generating-function approach.
By using a generating-function approach we are able to extend the applicability of the least-squares procedure originally presented in [5] to any optical system that can be described by a smooth generating function.