## Abstract

We present a unified mathematical framework for sixteen fundamental optical systems. The systems have a parallel or point source and a parallel, point, near-field or far-field target. These choices give eight configurations if we use reflectors only and take the minimum number of freeform surfaces required. Similarly, we get eight lens systems if we only use lens surfaces. The mathematical model for each system is based on Hamilton’s characteristic functions and conservation of luminous flux. Some configurations lead to standard or generalized Monge-Ampère equations. The remaining systems are described by so-called generated Jacobian equations.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

In this paper we present a unified mathematical framework for sixteen systems. These basic configurations have a parallel or point source and a parallel, point, near-field or far-field target. These choices give eight systems, if we use reflectors only and take the minimum number of freeform surfaces required. The systems are shown in Fig. 3. Similarly, we get eight basic lens systems if we only use lens surfaces, see Fig. 4. Our goal is to find the shape and location of the surfaces that convert a given source distribution into a desired target distribution. All systems can be modeled by a nonlinear partial differential equation (PDE). Many papers in freeform illumination optics cover only one particular optical system. Here we derive models for all sixteen using our unified framework.

We give a short overview of literature on finding freeform optical surfaces—a more detailed survey can be found in [1, Ch. 5]. First we consider papers with PDE-based models where a PDE for the surface is solved directly using finite difference discretizations. This approach is followed by Benamou, Froese and Oberman [2–9], Brix et al. [10,11], Caboussat et al. [12,13], Dean and Glowinski [14–16], Prins et al. [17–19], Romijn et al. [1,20–23], Wu et al. [24–30], and Yadav et al. [31–33]. Lakkis et al. [34–36] introduce a Galerkin-type finite element method that can solve PDEs of Monge-Ampère type.

A second approach is to solve the Monge-Kantorovich transport problem that corresponds to the design problem using optimization strategies, see Angenent et al. [37], Benamou et al. [38], and Haber et al. [39]. Doskolovich et al. reduce the transport problem to a linear assignment problem [40–42]. Oliker et al. have introduced the methods of supporting paraboloids and ellipsoids, now known as the supporting quadric method (SQM) [43–51]. SQM combines geometrical techniques with methods from calculus of variations. Instead of dealing directly with the complicated nonlinear PDEs, it produces weak solutions that are the envelope of a set of e.g. paraboloids of revolution. By construction there exists a supporting paraboloid at every point of the computed optical surface. The surfaces found are not necessarily differentiable. The method provides both proofs of existence of solutions as well as a way to compute them numerically.

Third are ray-mapping techniques. In these methods a ray mapping is computed which is subsequently used to compute the surface using the law of reflection or Snell’s law. This method is used by Bösel and Gross [52,53] and Feng et al. [54–56].

We have organized this paper as follows. In Section 2 we formulate the freeform design problem for a general optical system, introduce notation and recap the main properties of Hamilton’s characteristic functions. Next we present three mathematical models. In Sections 3 and 4 we discuss the eight reflector and eight lens systems, respectively. We sketch our least-squares algorithm to numerically solve the mathematical models in Section 5. We draw conclusions in Section 6.

## 2. Mathematical models

#### 2.1 Problem formulation

We consider a general optical system such as the one sketched in Fig. 1. The goal of freeform design is to find the shape and location of one or more surfaces that convert a given source distribution into a desired target distribution. The surfaces can be either reflector or lens surfaces.

The source distribution (emittance or intensity) is denoted by $f(\mathbf {x})$, where $\mathbf {x} \in X$ are either spatial or directional coordinates and $X$ is the source domain. Likewise, the target distribution (illuminance or intensity) is $g(\mathbf {y})$, $\mathbf {y} \in Y$, where $\mathbf {y} \in Y$ are spatial or directional coordinates and $Y$ is the target domain. The optical map $\mathbf {m} : X \to Y$ connects $\mathbf {x}$ with $\mathbf {y}$, so $\mathbf {y} = \mathbf {m}(\mathbf {x})$.

We use spatial coordinates for a parallel source, parallel target and near field target. We assume standard Cartesian coordinates in these cases. We use directional coordinates for a point source, point target and far field target. It is convenient to use so-called stereographic coordinates. These are a parametrization of the unit sphere: If we denote the direction vector of an emitted ray by $\hat {\mathbf {s}} = (s_1, s_2, s_3)$, then its stereographic projection onto the equatorial plane $z=0$ from the south pole is given by $\mathbf {x} = (x_1, x_2) = (s_1, s_2) / (1 + s_3)$, see Fig. 2. Note that the hat ($\hat {~}$) indicates that $\hat {\mathbf {s}}$ has length one. The inverse projection is given by $\hat {\mathbf {s}} = (2x_1, 2x_2, 1-|\mathbf {x}|^{2}) / (1 + |\mathbf {x}|^{2})$. In the same way, the direction vector $\hat {\mathbf {t}}=(t_1,t_2,t_3)$ of a target ray has stereographic projection $\mathbf {y}=(y_1,y_2)$. At the target, the projection can be either from the south or the north pole.

In Sections 2.3–2.5 we present models for the optical map. These models follow from the evaluation of one of the four Hamilton’s characteristic functions [57,58]. This gives a geometrical description of the optical system. We briefly state the definitions and main properties of the characteristic functions. Consider a light ray that intersects the source plane $z=z_s$ at (two-dimensional) spatial coordinates $\mathbf {q}_s$ and the target plane $z=z_t$ at spatial coordinates $\mathbf {q}_t$. The subscripts $s$ and $t$ refer to source and target. Hamilton’s point characteristic $V = V(\mathbf {q}_s, \mathbf {q}_t)$ is defined as the optical path length between $(\mathbf {q}_s,z_s)$ and $(\mathbf {q}_t,z_t)$, so $V = \int _C n \, \textrm {d} s$, with $C$ the path connecting the two points and $n$ the local refractive index. Let $\mathbf {p}_s$ be the vector with the first two components of the momentum at the source plane and, in the same way, $\mathbf {p}_t$ at the target plane. We have $\mathbf {p}_s = n (s_1, s_2)$ and $\mathbf {p}_t = n (t_1, t_2)$. The point characteristic $V$ satisfies

In our models, we use $V$ if both $\mathbf {x}$ and $\mathbf {y}$ are Cartesian coordinates, so if $\mathbf {x} = \mathbf {q}_s$ and $\mathbf {y} = \mathbf {q}_t$. The mixed characteristic $W$ is defined as $W(\mathbf {q}_s, \mathbf {p}_t) = V(\mathbf {q}_s, \mathbf {q}_t) - \mathbf {q}_t \boldsymbol {\cdot } \mathbf {p}_t$. Note that $W$ depends on two variables only: Using Eq. (1), we have ${\partial W}/{\partial \mathbf {q}_t} = {\partial V}/{\partial \mathbf {q}_t} - \mathbf {p}_t = \mathbf {0}.$ It satisfies

We will use $W$ if $\mathbf {x}$ is Cartesian ($\mathbf {x} = \mathbf {q}_s$) and $\mathbf {y}$ is the stereographic projection of $\hat {\mathbf {t}}$. The mixed characteristic of the second kind $W^{*}$ is defined as $W^{*}(\mathbf {p}_s, \mathbf {q}_t) = V(\mathbf {q}_s, \mathbf {q}_t) + \mathbf {q}_s \boldsymbol {\cdot } \mathbf {p}_s$. It depends on two variables only, because ${\partial W^{*}}/{\partial \mathbf {q}_s} = \mathbf {0}.$ It satisfies

We will use $W^{*}$ if $\mathbf {x}$ is the stereographic projection of $\hat {\mathbf {s}}$ and $\mathbf {y}$ is Cartesian ($\mathbf {y} = \mathbf {q}_t$). Finally, the angular characteristic $T$ is defined as $T(\mathbf {p}_s, \mathbf {p}_t) = V(\mathbf {q}_s, \mathbf {q}_t) + \mathbf {q}_s \boldsymbol {\cdot } \mathbf {p}_s - \mathbf {q}_t \boldsymbol {\cdot } \mathbf {p}_t$. It depends on two variables only, because ${\partial T}/{\partial \mathbf {q}_s} = {\partial T}/{\partial \mathbf {q}_t} = \mathbf {0}.$ It satisfies

We will use $T$ if $\mathbf {x}$ and $\mathbf {y}$ are the stereographic projections of $\hat {\mathbf {s}}$ and $\hat {\mathbf {t}}$, respectively.

#### 2.2 Conservation of luminous flux

Let $A$ be a part of the source domain, so $A \subset X$. Because luminous flux is conserved, we know that all light emitted from $A$ should reach the image set $\mathbf {m}(A)$ in the target domain. This leads to the conservation law

Here, $J_X(\mathbf {x})$ equals one if we use Cartesian coordinates in the source domain. If we use stereographic coordinates, it is the Jacobian of the stereographic projection of the direction vector $\hat {\mathbf {s}}$ of the emitted rays, so

In the same way $J_Y(\mathbf {y})$ equals one for Cartesian coordinates in the target domain and $J_Y(\mathbf {y}) = {4}/{ (1 + |\mathbf {y}|^{2})^{2} }$ for stereographic coordinates. $\textrm {D} \mathbf {m}$ is the Jacobian matrix of $\mathbf {m}$. We assume that $\det (\textrm {D} \mathbf {m}(\mathbf {x})) > 0$ in Eq. (5). For far-field targets we use the far-field approximation.

Because Eq. (5) needs to hold for any subset $A$ of $X$, we find the flux balance

Equation (7) will lead to elliptic PDEs. For all models, we impose the condition $\mathbf {m}(X) = Y$ stating that all light from the source domain must be transferred to the target domain.

#### 2.3 Model 1: standard Monge-Ampère equation

The first model we consider for the optical map will allow us to write $\mathbf {m}$ as the gradient of a function $u_1$, so $\mathbf {m} = \nabla u_1$. The function $u_1$ is related to the location of the optical surface. We will show that $u_1$ satisfies a second-order nonlinear PDE, the so-called Monge-Ampère equation.

A typical example of an optical system that can be described by this model is the parallel-to-far-field reflector sketched in Fig. 3(d). A detailed derivation of the model for this system can be found in [1, Sec. 3.2].

For the source we use Cartesian coordinates, so $J_X = 1$. For the target we use stereographic coordinates, so $J_Y(\mathbf {y}) = {4}/{ (1 + |\mathbf {y}|^{2})^{2} }$. We assume that the source domain is located at $z=0$ and the target domain at $z = -L$. Because emitted rays are parallel and perpendicular to $z=0$, we have $\mathbf {p}_s = \mathbf {0}$. With Eq. (2) it follows that Hamilton’s mixed characteristic $W$ is a function of $\mathbf {p}_t$ (and hence $\mathbf {y}$) only. We define the location of the reflector as $z = u_1(\mathbf {x})$. Evaluating $W$ gives $W(\mathbf {p}_t) = (1-t_3) u_1(\mathbf {x}) - \mathbf {x} \boldsymbol {\cdot } \mathbf {p}_t - L t_3$, which can be rewritten as

The first term in Eq. (8) depends on $\mathbf {y}$ only and we call it $u_2(\mathbf {y})$. For the right-hand side of Eq. (8), we use $\mathbf {p}_t/(1-t_3) = \mathbf {y}$. This leads to

A possible solution of Eq. (9) readsThe second equation in Eq. (10) implies $\nabla _{\mathbf {x}} \left ( -\mathbf {x} \boldsymbol {\cdot } \mathbf {y} + u_1(\mathbf {x}) \right ) = \mathbf {0}$, which leads to $\mathbf {y} = \mathbf {m}(\mathbf {x}) = \nabla u_1(\mathbf {x})$. Note that

Substitution of Eq. (11) into Eq. (7) gives a nonlinear elliptic PDE for $u_1$, viz.

We call this PDE the standard Monge-Ampère (SMA) equation.

#### 2.4 Model 2: generalized Monge-Ampère equation

In the second model, the optical map $\mathbf {m}$ is no longer the gradient of a function $u_1$. It is a (complicated) function of $\mathbf {x}$ and $\nabla u_1$, where $u_1$ is again related to the location of the optical surface. We will show that $u_1$ satisfies a generalized Monge-Ampère (GMA) equation.

A typical example of an optical system that can be described by this model is the point-to-far-field reflector sketched in Fig. 3(h). A detailed derivation of the model for this system can be found in [20].

For both the source and the target we use stereographic coordinates, so $J_X(\mathbf {x}) = {4}/{ (1 + |\mathbf {x}|^{2})^{2} }$ and $J_Y(\mathbf {y}) = {4}/{ (1 + |\mathbf {y}|^{2})^{2} }$. We assume that the source domain is located at $z=0$ and the target domain at $z = -L$. The reflector surface has parametrization $\mathbf {r}(\hat {\mathbf {s}}) = u(\hat {\mathbf {s}}) \hat {\mathbf {s}}$. Because we have a point source, we find $\mathbf {q}_s = \mathbf {0}$. With Eq. (4) it follows that Hamilton’s angular characteristic $T$ is a function of $\mathbf {p}_t$ (and hence $\mathbf {y}$) only. Evaluating $T$ gives $T(\mathbf {p}_t) = (1-\hat {\mathbf {s}}\boldsymbol {\cdot }\hat {\mathbf {t}}) u(\hat {\mathbf {s}}) - L t_3$, which can be rewritten as $\log (T(\mathbf {p}_t)+L t_3) = \log (1-\hat {\mathbf {s}}\boldsymbol {\cdot }\hat {\mathbf {t}}) + \log (u(\hat {\mathbf {s}}))$. Transforming to stereographic coordinates leads to

where $u_1(\mathbf {x}) = -\log (u(\hat {\mathbf {s}})) + \log ( 1 + |\mathbf {x}|^{2} )$, $u_2(\mathbf {y}) = -\log (T(\mathbf {p}_t)+L t_3) -\log (\frac 12( 1 + |\mathbf {y}|^{2} ))$ and $c(\mathbf {x}, \mathbf {y}) = -\log ( 1-2\mathbf {x}\boldsymbol {\cdot }\mathbf {y} + |\mathbf {x}|^{2} |\mathbf {y}|^{2} )$. The function $c$ is known as the cost function in optimal transport theory. Note that Eq. (13) is a generalization of Eq. (9). A possible solution of Eq. (13) readsThe second equation in Eq. (14) implies

In principle one can compute $\mathbf {m}$ from Eq. (15), but this is difficult and we proceed in a different manner. Substituting $\mathbf {y} = \mathbf {m}(\mathbf {x})$ and once more differentiating with respect to $\mathbf {x}$ leads to the following equation for $\textrm {D} \mathbf {m}$

Here $\textrm {D}^{2} u_1$ and $\textrm {D}_{\mathbf {x}\mathbf {x}} c$ are Hessian matrices, $\textrm {D}_{\mathbf {x}\mathbf {y}} c$ is the matrix with mixed second-order derivatives. We have now found the matrix equation

It follows that $\det (\mathbf {P}) = \det (\mathbf {C}) \det (\textrm {D} \mathbf {m})$. Note that $\mathbf {P}$ is a $2 \times 2$-matrix, so $\det (-\mathbf {P}) = \det (\mathbf {P})$. Using the definition of $\mathbf {P}$ and Eq. (7) we obtain the flux balance

#### 2.5 Model 3: generated Jacobian equation

In the third model, we cannot separate the variables as we did in the left-hand side of Eqs. (9) and (13). Instead we will find a more complicated relation between $u_1(\mathbf {x})$ and $u_2(\mathbf {y})$. We will show that $u_1$ satisfies a generated Jacobian equation (GJE).

A typical example of an optical system that can be described by this model is the parallel-to-near-field reflector sketched in Fig. 3(c). A detailed derivation of the model for this system can be found in [1, Sec. 3.3].

For both the source and target we use Cartesian coordinates, so $J_X = J_Y = 1$. We assume that the source domain is located at $z=0$ and the target domain at $z = -L$. Because emitted rays are parallel and perpendicular to $z=0$, we have $\mathbf {p}_s = \mathbf {0}$. With Eq. (1) it follows that Hamilton’s point characteristic $V$ is a function of $\mathbf {q}_t$ (and hence $\mathbf {y}$) only. We define the location of the reflector as $z = u_1(\mathbf {x})$. Evaluating $V$ gives $V(\mathbf {y}) = u_1(\mathbf {x}) + \sqrt { |\mathbf {x}-\mathbf {y}|^{2} + u_1(\mathbf {x})^{2} }$. Solving this equation for $u_1$ leads to the following generalization of Eq. (13):

where $u_2(\mathbf {y}) = V(\mathbf {y})$ and $G(\mathbf {x}, \mathbf {y}, v) = \frac 12 v - |\mathbf {x}-\mathbf {y}|^{2}/(2v)$. The function $G$ is the so-called generating function. Because $\partial G/\partial v > 0$, the function $G(\mathbf {x}, \mathbf {y},\cdot )$ has an inverse for fixed $\mathbf {x}$ and $\mathbf {y}$. We denote this inverse by $H(\mathbf {x}, \mathbf {y},\cdot )$, so $u_2(\mathbf {y}) = H(\mathbf {x}, \mathbf {y}, u_1(\mathbf {x}))$. A possible solution of Eq. (19) readsThe second equation in Eq. (20) implies

In principle one can compute $\mathbf {m}$ from Eq. (21), but this is difficult and we proceed in a different manner. We introduce $\widetilde H(\mathbf {x}, \mathbf {y}) := H(\mathbf {x}, \mathbf {y}, u_1(\mathbf {x}))$, substitute $\mathbf {y} = \mathbf {m}(\mathbf {x})$ and differentiate once more with respect to $\mathbf {x}$ to find the following equation for $\textrm {D} \mathbf {m}$

We have now found the matrix equation

It follows that $\det (\mathbf {P}) = \det (\mathbf {C}) \det (\textrm {D} \mathbf {m})$. Using the definition of $\mathbf {P}$ and Eq. (7) we obtain

## 3. Basic reflector systems

In this section we discuss the eight basic reflector systems sketched in Fig. 3. For all of the systems we assume that the source is located at $z=0$.

#### 3.1 Parallel-to-parallel reflector

We consider the parallel source to parallel target reflector system sketched in Fig. 3(a). This system can be described by Model 1 from Section 2.3. A detailed derivation of the model for this system can be found in [32].

For both the source and target we use Cartesian coordinates, so $J_X = J_Y = 1$. We assume that the target is located at $z=L$. The first reflector has equation $z=u_1(\mathbf {x})$. The second reflector has equation $z=L+u_2(\mathbf {y})$. Because both emitted and target rays are parallel and perpendicular to $z=0$ and $z=L$, we have $\mathbf {p}_s = \mathbf {p}_t = \mathbf {0}$. With Eq. (1) it follows that Hamilton’s point characteristic $V$ is constant. This is also known as the principle of Malus and Dupin [57, Sec. 3.3.3]. If we introduce the reduced optical path length $\beta =V-L$, we get Eq. (13) with $c(\mathbf {x}, \mathbf {y}) := |\mathbf {x}-\mathbf {y}|^{2}/(2\beta ) - \frac 12\beta - L$. Equation (15) leads to $(\mathbf {x}-\mathbf {y})/\beta + \nabla u_1(\mathbf {x}) = \mathbf {0}$. We find

Substitution of Eq. (25) into Eq. (7) gives the SMA equation

#### 3.2 Parallel-to-point reflector

In this section we consider the parallel source to point target reflector system sketched in Fig. 3(b). This system can be described by Model 2 from Section 2.4. Note that the system is the same as the point source to parallel target system, that we will consider in Section 3.5, when we interchange source and target coordinates. A detailed derivation of the model for the latter system can be found in [59].

For the source we use Cartesian coordinates, so $J_X = 1$. For the target we use stereographic coordinates, so $J_Y(\mathbf {y}) = {4}/{ (1 + |\mathbf {y}|^{2})^{2} }$. We assume that the target is located at $z=L$. Because emitted rays are parallel and perpendicular to $z=0$, we have $\mathbf {p}_s = \mathbf {0}$. The target is a point, so $\mathbf {q}_t = \mathbf {0}$. With Eq. (2) it follows that Hamilton’s mixed characteristic $W$ is constant. If we introduce the reduced optical path length $\beta =W-L$, we get Eq. (13) with

Repeating the steps described in Section 2.4 leads to the GMA equation

#### 3.3 Parallel-to-near-field reflector

We consider the parallel source to near field reflector system sketched in Fig. 3(c). This system can be described by Model 3 from Section 2.5, and in fact it is the system discussed in that section.

#### 3.4 Parallel-to-far-field reflector

We consider the parallel source to far field reflector system sketched in Fig. 3(d). This system can be described by Model 1 from Section 2.3, and in fact it is the system discussed in that section.

#### 3.5 Point-to-parallel reflector

In this section we consider the point source to parallel target reflector system sketched in Fig. 3(e). This system can be described by Model 2 from Section 2.4. A detailed derivation of the model for this system can be found in [59]. Note that the system is the same as the parallel source to point target system, described in Section 3.2, when we interchange source and target coordinates.

For the source we use stereographic coordinates, so $J_X(\mathbf {x}) = {4}/{ (1 + |\mathbf {x}|^{2})^{2} }$. For the target we use Cartesian coordinates, so $J_Y = 1$. We assume that the target is located at $z=L$. Because we have a point source, we find $\mathbf {q}_s = \mathbf {0}$. Target rays are parallel and perpendicular to $z=L$, so $\mathbf {p}_t = \mathbf {0}$. With Eq. (3) it follows that Hamilton’s mixed characteristic $W^{*}$ is constant. If we introduce the reduced optical path length $\beta =W^{*}-L$, we get Eq. (13) with

Repeating the steps in Section 2.4 leads to the GMA equation

#### 3.6 Point-to-point reflector

In this section we consider the point source to point target reflector system sketched in Fig. 3(f). This system can be described by Model 2 from Section 2.4. A detailed derivation of the model for this system can be found in [33, Sec. 3.4].

For both the source and the target we use stereographic coordinates, so $J_X(\mathbf {x}) = {4}/{ (1 + |\mathbf {x}|^{2})^{2} }$ and $J_Y(\mathbf {y}) = {4}/{ (1 + |\mathbf {y}|^{2})^{2} }$. We assume that the target is located at $z=L$. Because we have both a point source and a point target, we find $\mathbf {q}_s = \mathbf {q}_t = \mathbf {0}$. With Eq. (4) it follows that Hamilton’s angular characteristic $T$ is constant. We get Eq. (13) with

Repeating the steps in Section 2.4 leads to the GMA equation

#### 3.7 Point-to-near-field reflector

In this section we consider the point source to near field reflector system sketched in Fig. 3(g). This system can be described by Model 3 from Section 2.5.

For the source we use stereographic coordinates, so $J_X(\mathbf {x}) = {4}/{ (1 + |\mathbf {x}|^{2})^{2} }$. For the target we use Cartesian coordinates, so $J_Y = 1$. We assume that the target is located at $z=L$. Because we have a point source, we find $\mathbf {q}_s = \mathbf {0}$. With Eq. (3) it follows that Hamilton’s mixed characteristic $W^{*}$ is a function of $\mathbf {q}_t$ (and hence $\mathbf {y}$) only. Evaluating $W^{*}$ we obtain Eq. (19) with

Repeating the steps in Section 2.5 we find the GJE

#### 3.8 Point-to-far-field reflector

In this section we consider the point source to far field reflector system sketched in Fig. 3(h). This system can be described by Model 2 from Section 2.4, and in fact it is the system discussed in that section.

## 4. Basic lens systems

In this section we discuss the eight basic lens systems sketched in Fig. 4. For all of the systems we assume that the source is located at $z=0$, the target at $z=L$ and that the lens has refractive index $n > 1$.

#### 4.1 Parallel-to-parallel lens

We consider the parallel source to parallel target lens system sketched in Fig. 4(a). This system can be described by Model 2 from Section 2.4. A detailed derivation of the model for this system can be found in [31].

For both the source and target we use Cartesian coordinates, so $J_X = J_Y = 1$. Because both emitted and target rays are parallel and perpendicular to $z=0$ and $z=L$, we have $\mathbf {p}_s = \mathbf {p}_t = \mathbf {0}$. With Eq. (1) it follows that Hamilton’s point characteristic $V$ is constant. If we introduce the reduced optical path length $\beta =V-L$, we get Eq. (13) with

Repeating the steps in Section 2.4 leads to the GMA equation

#### 4.2 Parallel-to-point lens

In this section we consider the parallel source to point target lens system sketched in Fig. 4(b). This system can be described by Model 3 from Section 2.5. Note that the system is the same as the point source to parallel target system, that we will consider in Section 4.5, when we interchange source and target coordinates.

For the source we use Cartesian coordinates, so $J_X = 1$. For the target we use stereographic coordinates, so $J_Y(\mathbf {y}) = {4}/{ (1 + |\mathbf {y}|^{2})^{2} }$. Because emitted rays are parallel and perpendicular to $z=0$, we have $\mathbf {p}_s = \mathbf {0}$. The target is a point, so $\mathbf {q}_t = \mathbf {0}$. With Eq. (2) it follows that Hamilton’s mixed characteristic $W$ is constant. If we introduce the reduced optical path length $\beta =W-L$, we get Eq. (19) with

#### 4.3 Parallel-to-near-field lens

We consider the parallel source to near field lens system sketched in Fig. 4(c). Note that the first lens surface is flat and perpendicular to the emitted rays, so it does not change their direction. This system can be described by Model 3 from Section 2.5.

For both the source and target we use Cartesian coordinates, so $J_X = J_Y = 1$. Because emitted rays are parallel and perpendicular to $z=0$, we have $\mathbf {p}_s = \mathbf {0}$. With Eq. (1) it follows that Hamilton’s point characteristic $V$ is a function of $\mathbf {q}_t$ (and hence $\mathbf {y}$) only. Evaluating $V$ we obtain Eq. (19) with

Repeating the steps in Section 2.5 we find the GJE

#### 4.4 Parallel-to-far-field lens

We consider the parallel source to far field lens system sketched in Fig. 4(d). Note that the first lens surface is flat and perpendicular to the emitted rays, so it does not change their direction. This system can be described by Model 1 from Section 2.3. A detailed derivation of the model for this system can be found in [33, Sec. 3.1.2].

For the source we use Cartesian coordinates, so $J_X = 1$. For the target we use stereographic coordinates, so $J_Y(\mathbf {y}) = {4}/{ (1 + |\mathbf {y}|^{2})^{2} }$. Because emitted rays are parallel and perpendicular to $z=0$, we have $\mathbf {p}_s = \mathbf {0}$. With Eq. (2) it follows that Hamilton’s mixed characteristic $W$ is a function of $\mathbf {p}_t$ (and hence $\mathbf {y}$) only. Evaluating $W$ we get Eq. (9) where $u_1(\mathbf {x})$ defines the location of the freeform lens surface via $z = u_1(\mathbf {x})$. Repeating the steps in Section 2.3 we find the SMA Eq. (12).

#### 4.5 Point-to-parallel lens

In this section we consider the point source to parallel target lens system sketched in Fig. 4(e). This system can be described by Model 3 from Section 2.5. Note that the system is the same as the parallel source to point target system, described in Section 4.2, when we interchange source and target coordinates.

For the source we use stereographic coordinates, so $J_X(\mathbf {x}) = {4}/{ (1 + |\mathbf {x}|^{2})^{2} }$. For the target we use Cartesian coordinates, so $J_Y = 1$. Because we have a point source, we find $\mathbf {q}_s = \mathbf {0}$. Target rays are parallel and perpendicular to $z=L$, so $\mathbf {p}_t = \mathbf {0}$. With Eq. (3) it follows that Hamilton’s mixed characteristic $W^{*}$ is constant. If we introduce the reduced optical path length $\beta =W^{*}-L$, we get Eq. (19) with

#### 4.6 Point-to-point lens

In this section we consider the point source to point target lens system sketched in Fig. 4(f). This system can be described by Model 3 from Section 2.5.

For both the source and the target we use stereographic coordinates, so $J_X(\mathbf {x}) = {4}/{ (1 + |\mathbf {x}|^{2})^{2} }$ and $J_Y(\mathbf {y}) = {4}/{ (1 + |\mathbf {y}|^{2})^{2} }$. Because we have both a point source and a point target, we find $\mathbf {q}_s = \mathbf {q}_t = \mathbf {0}$. With Eq. (4) it follows that Hamilton’s angular characteristic $T$ is constant. We get Eq. (19) with

#### 4.7 Point-to-near-field lens

In this section we consider the point source to near field lens system sketched in Fig. 4(g). This system can be described by Model 3 from Section 2.5.

For the source we use stereographic coordinates, so $J_X(\mathbf {x}) = {4}/{ (1 + |\mathbf {x}|^{2})^{2} }$. For the target we use Cartesian coordinates, so $J_Y = 1$. The first lens surface is spherical and does not change the direction of the rays. The second lens surface is freeform. Because we have a point source, we find $\mathbf {q}_s = \mathbf {0}$. With Eq. (3) it follows that Hamilton’s mixed characteristic $W^{*}$ is a function of $\mathbf {q}_t$ (and hence $\mathbf {y}$) only. Evaluating $W^{*}$ we obtain Eq. (19) with

#### 4.8 Point-to-far-field lens

In this section we consider the point source to far field lens system sketched in Fig. 4(h). This system can be described by Model 2 from Section 2.4. A detailed derivation of the model for this system can be found in [21].

For both the source and the target we use stereographic coordinates, so $J_X(\mathbf {x}) = {4}/{ (1 + |\mathbf {x}|^{2})^{2} }$ and $J_Y(\mathbf {y}) = {4}/{ (1 + |\mathbf {y}|^{2})^{2} }$. The first lens surface is spherical and does not change the direction of the rays. The second lens surface is freeform and has parametrization $\mathbf {r}(\hat {\mathbf {s}}) = u(\hat {\mathbf {s}}) \hat {\mathbf {s}}$. Because we have a point source, we find $\mathbf {q}_s = \mathbf {0}$. With Eq. (4) it follows that Hamilton’s angular characteristic $T$ is a function of $\mathbf {p}_t$ (and hence $\mathbf {y}$) only. Evaluating $T$ we get Eq. (13) with

Repeating the steps in Section 2.4 leads to the GMA Eq. (18).

## 5. Numerical method for finding the optical surface(s)

Our numerical method for finding the optical surface(s) is an iterative least-squares solver. It is based on the PDE of SMA, GMA or GJE type together with the condition $\mathbf {m}(X) = Y$. Our implementation is based on the following breakdown in substeps. We compute $\mathbf {P}$, $\mathbf {b}$ and $\mathbf {m}$ such that

Here, $\mathbf {b}$ maps the boundary of the source domain to the boundary of the target domain. Equations (48) are solved in an iterative procedure: We let $n = 0$ and choose an initial guess $\mathbf {m}^{0}$. Next, let $J_I(\mathbf {m}, \mathbf {P}) = \tfrac 12 \iint _X \| \mathbf {C} \textrm {D} \mathbf {m} - \mathbf {P} \|^{2} \, \textrm {d}\mathbf {x}$ and compute $\mathbf {P}^{n+1}$ such that it minimizes $J_I(\mathbf {m}^{n}, \mathbf {P})$. This leads to a constrained minimization problem—the contraint is Eq. (48b). To enforce the boundary condition Eq. (48c), we define $J_B(\mathbf {m}, \mathbf {b}) = \tfrac 12 \int _{\partial X} | \mathbf {m} - \mathbf {b} |^{2} \, \textrm {d} s$ and compute $\mathbf {b}^{n+1}$ such that it minimizes $J_B(\mathbf {m}^{n}, \mathbf {b})$. We finally introduce $J(\mathbf {m}, \mathbf {P}, \mathbf {b}) = \alpha \, J_I(\mathbf {m}, \mathbf {P}) + (1-\alpha ) \, J_B(\mathbf {m}, \mathbf {b})$ where $0<\alpha <1$ and compute $\mathbf {m}^{n+1}$ such that it minimizes $J(\mathbf {m}, \mathbf {P}^{n+1}, \mathbf {b}^{n+1})$. This leads to an elliptic PDE for $\mathbf {m}$ that we solve with a finite volume method.

For Models 1 and 2, we repeat the steps above until convergence. Upon convergence of the iterative procedure, the location of the surface is calculated from the optical map using Eq. (15). For Model 3, we need to do an additional step during the iteration to compute $u_1$. To do so, we use Eq. (21) and let $I(u_1, \mathbf {m}) = \frac 12 \iint _X \left | \nabla _{\mathbf {x}} H(\mathbf {x}, \mathbf {m}, u_1) + H_v(\mathbf {x}, \mathbf {m}, u_1) \nabla u_1 \right |^{2} \, \textrm {d}\mathbf {x}$, where $H_v$ is the partial derivative of $H$ with respect to its third argument. We compute $u_1^{n+1}$ such that it minimizes $I(u_1, \mathbf {m}^{n+1})$.

For details on the numerical algorithm for SMA and GMA equations, see [20]. A full description on the generalization to GJEs can be found in [23]. Numerical results computed with this solver can be found in [18–23,31,32,59].

## 6. Conclusions

In this paper we presented a unified mathematical framework for sixteen basic optical systems, eight reflector and eight lens systems. Each system has a parallel or point source and a parallel, point, near-field or far-field target. The mathematical model for each system is based on Hamilton’s characteristic functions and conservation of luminous flux.

Three systems lead to the standard Monge-Ampère (SMA) equation (Model 1). This holds true for the parallel-to-parallel reflector, the parallel-to-far-field reflector and the parallel-to-far-field lens system.

More complicated optical systems cannot be described by Model 1 and need a generalization. Six of them give a generalized Monge-Ampère (GMA) equation (Model 2). Systems in this category are the parallel-to-point reflector, the point-to-parallel reflector, the point-to-point reflector, the point-to-far-field reflector, the parallel-to-parallel lens and the point-to-far-field lens.

The remaining seven systems cannot be described by either Model 1 or 2. They need a yet more general formulation using so-called generating functions. This in turn leads to generated Jacobian equations (GJEs, Model 3). All systems with near-field targets belong to this category.

Note that the three models are proper generalizations in the sense that all Model 1 systems can also be described by Model 2, and Model 2 systems can be formulated using generating functions (Model 3) too.

## Disclosures

The authors declare no conflicts of interest.

## Data availability

No data were generated or analyzed in the presented research.

## References

**1. **L. B. Romijn, “Generated Jacobian equations in freeform optical design,” ’ Ph.D. thesis, Eindhoven University of Technology, Eindhoven, the Netherlands (2021).

**2. **J.-D. Benamou, B. D. Froese, and A. M. Oberman, “Numerical solution of the optimal transportation problem using the Monge-Ampère equation,” J. Comput. Phys. **260**, 107–126 (2014). [CrossRef]

**3. **B. D. Froese, “A numerical method for the elliptic Monge-Ampère equation with transport boundary conditions,” SIAM J. Sci. Comput. **34**(3), A1432–A1459 (2012). [CrossRef]

**4. **B. D. Froese and A. M. Oberman, “Convergent filtered schemes for the Monge–Ampère partial differential equation,” SIAM J. Numer. Anal. **51**(1), 423–444 (2013). [CrossRef]

**5. **A. M. Oberman, “Convergent difference schemes for degenerate elliptic and parabolic equations: Hamilton–Jacobi equations and free boundary problems,” SIAM J. Numer. Anal. **44**(2), 879–895 (2006). [CrossRef]

**6. **A. M. Oberman, “Wide stencil finite difference schemes for the elliptic Monge-Ampère equation and functions of the eigenvalues of the Hessian,” Discret. Contin. Dyn. Syst. B **10**(1), 221–238 (2008). [CrossRef]

**7. **B. D. Froese and A. M. Oberman, “Fast finite difference solvers for singular solutions of the elliptic Monge–Ampère equation,” J. Comput. Phys. **230**(3), 818–834 (2011). [CrossRef]

**8. **J.-D. Benamou, B. D. Froese, and A. M. Oberman, “Two numerical methods for the elliptic Monge-Ampère equation,” Math. Model. Numer. Anal. **44**(4), 737–758 (2010).

**9. **B. D. Froese and A. M. Oberman, “Convergent finite difference solvers for viscosity solutions of the elliptic Monge–Ampère equation in dimensions two and higher,” SIAM J. Numer. Anal. **49**(4), 1692–1714 (2011). [CrossRef]

**10. **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). [CrossRef]

**11. **K. Brix, Y. Hafizogullari, and A. Platen, “Solving the Monge-Ampère equations for the inverse reflector problem,” Math. Models Methods Appl. Sci. **25**(05), 803–837 (2015). [CrossRef]

**12. **A. Caboussat, R. Glowinski, and D. C. Sorensen, “A least-squares method for the numerical solution of the Dirichlet problem for the elliptic Monge-Ampère equation in dimension two,” ESAIM: Control. Optim. Calc. Var. **19**(3), 780–810 (2013). [CrossRef]

**13. **A. Caboussat, R. Glowinski, and D. Gourzoulidis, “A least-squares/relaxation method for the numerical solution of the three-dimensional elliptic Monge-Ampère equation,” J. Sci. Comput. **77**(1), 53–78 (2018). [CrossRef]

**14. **E. J. Dean and R. Glowinski, “Numerical solution of the two-dimensional elliptic Monge-Ampère equation with Dirichlet boundary conditions: An augmented Lagrangian approach,” Comptes rendus Mathématique **336**(9), 779–784 (2003). [CrossRef]

**15. **E. J. Dean and R. Glowinski, “Numerical solution of the two-dimensional elliptic Monge-Ampère equation with Dirichlet boundary conditions: A least-squares approach,” Comptes rendus Mathématique **339**(12), 887–892 (2004). [CrossRef]

**16. **E. J. Dean and R. Glowinski, “Numerical methods for fully nonlinear elliptic equations of the Monge-Ampère type,” Comput. Method. Appl. M. **195**(13-16), 1344–1386 (2006). [CrossRef]

**17. **C. R. Prins, “Inverse methods for illumination optics,” Ph.D. thesis, Eindhoven University of Technology, Eindhoven, the Netherlands (2014).

**18. **C. R. Prins, R. Beltman, J. H. M. ten Thije Boonkkamp, W. L. IJzerman, and T. W. Tukker, “A least-squares method for optimal transport using the Monge-Ampère equation,” SIAM J. Sci. Comput. **37**(6), B937–B961 (2015). [CrossRef]

**19. **C. R. Prins, J. H. M. ten Thije Boonkkamp, J. van Roosmalen, W. L. IJzerman, and T. W. Tukker, “A Monge-Ampère-solver for free-form reflector design,” SIAM J. on Sci. Comput. **36**(3), B640–B660 (2014). [CrossRef]

**20. **L. B. Romijn, J. H. M. ten Thije Boonkkamp, and W. L. IJzerman, “Inverse reflector design for a point source and far-field target,” J. Comput. Phys. **408**, 109283 (2020). [CrossRef]

**21. **L. B. Romijn, J. H. M. ten Thije Boonkkamp, and W. L. IJzerman, “Freeform lens design for a point source and far-field target,” J. Opt. Soc. Am. A **36**(11), 1926–1939 (2019). [CrossRef]

**22. **L. B. Romijn, J. H. M. ten Thije Boonkkamp, M. J. H. Anthonissen, and W. L. IJzerman, “An iterative least-squares method for generated Jacobian equations in freeform optical design,” SIAM J. Sci. Comput. **43**(2), B298–B322 (2021). [CrossRef]

**23. **L. B. Romijn, J. H. M. ten Thije Boonkkamp, M. J. H. Anthonissen, and W. L. IJzerman, “Generating-function approach for double freeform lens design,” J. Opt. Soc. Am. A **38**(3), 356–368 (2021). [CrossRef]

**24. **R. Wu, L. Xu, P. Liu, Y. Zhang, Z. Zheng, H. Li, and X. Liu, “Freeform illumination design: A nonlinear boundary problem for the elliptic Monge-Ampère equation,” Opt. Lett. **38**(2), 229–231 (2013). [CrossRef]

**25. **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). [CrossRef]

**26. **Y. Zhang, R. Wu, P. Liu, Z. Zheng, H. Li, and X. Liu, “Double freeform surfaces design for laser beam shaping with Monge-Ampère equation method,” Opt. Commun. **331**, 297–305 (2014). [CrossRef]

**27. **R. Wu and H. Hua, “Direct design of aspherical lenses for extended non-Lambertian sources in three-dimensional rotational geometry,” Opt. Express **24**(2), 1017–1030 (2016). [CrossRef]

**28. **R. Wu, C. Y. Huang, X. Zhu, H.-N. Cheng, and R. Liang, “Direct three-dimensional design of compact and ultra-efficient freeform lenses for extended light sources,” Optica **3**(8), 840–843 (2016). [CrossRef]

**29. **S. Chang, R. Wu, L. An, and Z. Zheng, “Design beam shapers with double freeform surfaces to form a desired wavefront with prescribed illumination pattern by solving a Monge-Ampère type equation,” J. Opt. **18**(12), 125602 (2016). [CrossRef]

**30. **L. Yang, Y. Liu, Z. Ding, J. Zhang, X. Tao, Z. Zheng, and R. Wu, “Design of freeform lenses for illuminating hard-to-reach areas through a light-guiding system,” Opt. Express **28**(25), 38155–38168 (2020). [CrossRef]

**31. **N. K. Yadav, J. H. M. ten Thije Boonkkamp, and W. L. IJzerman, “A Monge-Ampère problem with non-quadratic cost function to compute freeform lens surfaces,” J. Sci. Comput. **80**(1), 475–499 (2019). [CrossRef]

**32. **N. K. Yadav, L. B. Romijn, J. H. M. ten Thije Boonkkamp, and W. L. IJzerman, “A least-squares method for the design of two-reflector optical systems,” JPhys Photonics **1**(3), 034001 (2019). [CrossRef]

**33. **N. K. Yadav, “Monge-Ampère problems with non-quadratic cost function. Application to freeform optics,’ Ph.D. thesis, Eindhoven University of Technology, Eindhoven, the Netherlands (2018).

**34. **O. Lakkis and T. Pryer, “A finite element method for second order nonvariational elliptic problems,” SIAM J. Sci. Comput. **33**(2), 786–801 (2011). [CrossRef]

**35. **O. Lakkis and T. Pryer, “A finite element method for nonlinear elliptic problems,” SIAM J. Sci. Comput. **35**(4), A2025–A2045 (2013). [CrossRef]

**36. **E. L. Kawecki, O. Lakkis, and T. Pryer, “A finite element method for the Monge-Ampère equation with transport boundary conditions,” arXiv preprint arXiv:1807.03535 (2018).

**37. **S. B. Angenent, S. Haker, and A. R. Tannenbaum, “Minimizing flows for the Monge-Kantorovich problem,” SIAM J. Math. Anal. **35**(1), 61–97 (2003). [CrossRef]

**38. **J.-D. Benamou, G. Carlier, M. Cuturi, L. Nenna, and G. Peyré, “Iterative Bregman projections for regularized transportation problems,” SIAM J. Sci. Comput. **37**(2), A1111–A1138 (2015). [CrossRef]

**39. **E. Haber, T. Rehman, and A. R. Tannenbaum, “An efficient numerical method for the solution of the *L*_{2} optimal mass transfer problem,” SIAM J. Sci. Comput. **32**(1), 197–211 (2010). [CrossRef]

**40. **L. L. Doskolovich, D. A. Bykov, E. S. Andreev, E. A. Bezus, and V. I. Oliker, “Designing double freeform surfaces for collimated beam shaping with optimal mass transportation and linear assignment problems,” Opt. Express **26**(19), 24602–24613 (2018). [CrossRef]

**41. **L. L. Doskolovich, D. A. Bykov, A. A. Mingazov, and E. A. Bezus, “Optimal mass transportation and linear assignment problems in the design of freeform refractive optical elements generating far-field irradiance distributions,” Opt. Express **27**(9), 13083–13097 (2019). [CrossRef]

**42. **D. A. Bykov, L. L. Doskolovich, A. A. Mingazov, E. A. Bezus, and N. L. Kazanskiy, “Linear assignment problem in the design of freeform refractive optical elements generating prescribed irradiance distributions,” Opt. Express **26**(21), 27812–27825 (2018). [CrossRef]

**43. **V. I. Oliker, “On reconstructing a reflecting surface from the scattering data in the geometric optics approximation,” Inverse Probl. **5**(1), 51–65 (1989). [CrossRef]

**44. **S. A. Kochengin and V. I. Oliker, “Determination of reflector surfaces from near-field scattering data,” Inverse Probl. **13**(2), 363–373 (1997). [CrossRef]

**45. **S. A. Kochengin and V. I. Oliker, “Determination of reflector surfaces from near-field scattering data II. Numerical solution,” Numer. Math. **79**(4), 553–568 (1998). [CrossRef]

**46. **T. Glimm and V. I. Oliker, “Optical design of single reflector systems and the Monge-Kantorovich mass transfer problem,” J. Math. Sci. **117**(3), 4096–4108 (2003). [CrossRef]

**47. **S. A. Kochengin and V. I. Oliker, “Computational algorithms for constructing reflectors,” Comput. Vis. Sci. **6**(1), 15–21 (2003). [CrossRef]

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

**49. **V. I. Oliker, J. Rubinstein, and G. M. Wolansky, “Supporting quadric method in optical design of freeform lenses for illumination control of a collimated light,” Adv. Appl. Math. **62**, 160–183 (2015). [CrossRef]

**50. **V. I. Oliker, “Controlling light with freeform multifocal lens designed with supporting quadric method (SQM),” Opt. Express **25**(4), A58–A72 (2017). [CrossRef]

**51. **V. I. Oliker, L. L. Doskolovich, and D. A. Bykov, “Beam shaping with a plano-freeform lens pair,” Opt. Express **26**(15), 19406–19419 (2018). [CrossRef]

**52. **C. Bösel and H. Gross, “Ray mapping approach for the efficient design of continuous freeform surfaces,” Opt. Express **24**(13), 14271–14282 (2016). [CrossRef]

**53. **C. Bösel and H. Gross, “Double freeform illumination design for prescribed wavefronts and irradiances,” J. Opt. Soc. Am. A **35**(2), 236–243 (2018). [CrossRef]

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

**55. **Z. Feng, B. D. Froese, C.-Y. Huang, D. Ma, and R. Liang, “Creating unconventional geometric beams with large depth of field using double freeform-surface optics,” Appl. Opt. **54**(20), 6277–6281 (2015). [CrossRef]

**56. **Z. Feng, D. Cheng, and Y. Wang, “Iterative wavefront tailoring to simplify freeform optical design for prescribed irradiance,” Opt. Lett. **44**(9), 2274–2277 (2019). [CrossRef]

**57. **M. Born and E. Wolf, * Principles of Optics: 60th Anniversary Edition* (Cambridge University, 2019), 7th ed.

**58. **R. K. Luneburg and M. Herzberger, * Mathematical Theory of Optics* (University of California, 1964).

**59. **A. H. van Roosmalen, M. J. H. Anthonissen, W. L. IJzerman, and J. H. M. ten Thije Boonkkamp, “Design of a freeform two-reflector system to collimate and shape a point source distribution,” Opt. Express **29**(16), 25605–25625 (2021). [CrossRef]