## Abstract

In this work, a new two-dimensional optics design method is proposed that enables the coupling of three ray sets with two lens surfaces. The method is especially important for optical systems designed for wide field of view and with clearly separated optical surfaces. Fermat’s principle is used to deduce a set of functional differential equations fully describing the entire optical system. The presented general analytic solution makes it possible to calculate the lens profiles. Ray tracing results for calculated 15^{th} order Taylor polynomials describing the lens profiles demonstrate excellent imaging performance and the versatility of this new analytic design method.

© 2012 OSA

## 1. Introduction

In geometrical optics, the modern formulation of Fermat’s principle states that rays of light traverse the path of stationary optical length

between two points*A*and

*B*through media with refractive index distribution

*n*. As done in most textbooks on optics, it can be used to describe the properties of light rays refracted through different media, reflected off mirrors or undergoing total internal reflection. Indeed, all known laws of geometrical optics, lens design and aberrations will be consequences of the analytic properties of solutions to Fermat’s principle [1]. However, its mathematical complexity could give rise to the impression that it is impractical to be used for optics design. In fact, conventional optics design is based on minimizing a chosen merit function which quantifies the system’s performance for a defined sets of rays. In case of imaging applications, this can be merit functions such as the sum of the squares of certain aberrations or the RMS blur spot at the image plane [2]. In case of nonimaging optics, different merit functions such as the contrast ratio at a receiver plane or the optical collection efficiency could be chosen.

For many optical design problems, perfect solutions normally do not exist. Therefore, the design strategy begins with a parameterized description of the (refractive or reflective) optical surfaces. By using multi-parametric optimization (a common tool of any optical design software), these optimization algorithms normally start from an initial set of parameters and end with a final optimized set of parameters that minimizes a defined merit function. Since this is often a non-convex optimization problem, it cannot be guaranteed that these algorithms find the global minimum.

In contrast, direct optics design methods do not necessarily follow this design strategy. Usually, a key objective is to find a set of unrestricted free-form surfaces designed to have specific predefined characteristics. Where unrestricted means that the optical surfaces can have any shape to fulfill all imposed requirements. Such an approach can help to reduce the needed number of optical elements to a minimum or offer a way to design very compact optical systems. A basic example of such a design task is to focus light coming from a point in object space onto a point in the image space which can be achieved by using a Cartesian oval [3]. In order to focus light coming from an additional object point, one surface is no longer sufficient, two surfaces are needed. In general, an optical system consisting of *N* optical surfaces can couple *N* sets of rays for which specific conditions are imposed.

In case of optical systems designed for wide field of view and with at least one surface far from the aperture stop, it will be shown in this paper that it is possible to couple more than two ray sets with only two lens surfaces. However, this can only be achieved if different ray sets use different portions of the lens surfaces. Optical systems, where different incident directions use different portions of lens’ surfaces, are widely known. Field-flattener lenses are used in binocular designs and in astronomic telescopes to improve edge sharpness and lower distortion. Aperture stops in imaging systems often target the same objective. Based on a very basic example of a single thick lens, the Simultaneous Multiple Surfaces (SMS) design method will serve as a starting point in Sec. 2. It will help to provide a better understanding of the practical implications on the design process of an increased lens thickness and a wider field of view. This foundation will lead to an irreducible design representation in Sec. 3, based on the local coupling of two ray sets with two surface parts to achieve global coupling of three ray sets with two smooth surfaces. A first implemented numerical algorithm demonstrates the existence of solutions and it is used to introduce and establish the concept of convergence points which were first introduced to design aspheric V-groove reflectors [4]. Fermat’s principle is then used to deduce a system of functional differential equations in Sec. 4, fully describing the lens’ profiles. The transformation of these functional differential equations into an algebraic linear system of equations allows the successive calculation of the Taylor series coefficients up to an arbitrary order. Ray tracing analysis in Sec. 5 is then used to demonstrate the versatility of this analytic design approach by explicitly evaluating Taylor polynomials of 15^{th} order.

## 2. Initial degrees of freedom for SMS2D designs

The Simultaneous Multiple Surfaces (SMS) method has proven to be very versatile in various applications. SMS surfaces are piecewise curves made of several portions of generalized Cartesian ovals that map initial ray sets to final ray sets. It involves the simultaneous calculation of *N* optical surfaces using *N* one-parameter ray sets for which specific conditions connect the initial with the final ray sets [3]. A particular formulation of the SMS2D method for imaging systems comprises perfect imaging of *N* ray sets at the correspondent *N* image points. The SMS2D method offers the flexibility to choose the ray sets and their associated image points [5]. SMS optics are calculated by applying a constant optical path length for each coupled ray set. For two design angles of opposite sign the overall symmetry implies an identical optical path length for both ray sets. The optical path length can be determined by choosing one initial point on each lens profile. Figure 1(a)–(c) show different SMS2D designs for *θ* = ±5° design angle and increasing lens thicknesses from left to right. Similarly, Fig. 1(d)–(f) shows designs with increasing lens thicknesses for *θ* = ±10° design angle.

The bold light rays highlight the initial points used to determine the optical path length. Additional rays alternately construct the final lens profiles by applying the constant optical path length condition. This design process only provides a set of discrete points and correspondent normal vectors. An *SMS skinning* process can finally be used to fill the gaps between the separated points, fully defining both optical surfaces. A segment can be interpolated between the two symmetric adjacent points of the bottom surface, as indicated as solid lens profile in Fig. 1. As a boundary condition, the segment has to match point coordinates and normal vectors at the initial points, which requires at least a 2^{nd} order polynomial function. However, in principle, any analytical function satisfying this boundary condition can be applied - which correspondents to an infinite parameter space for the initial segment problem. For small design angles and moderate lens thicknesses, the initial segment represents only a small fraction of the entire lens profiles. However, with increasing design angles and/or distances between the optical surfaces the relative importance of the initial segment increases as well. In extreme cases, the initial segment can completely define the full lens profiles as shown in Fig. 1(f). A selected 2^{nd} order polynomial segment may guarantee coupling of two ray sets with chosen design angles of opposite sign, but it does not make any further use of the full potential offered by an unrestricted initial segment satisfying the boundary condition. Therefore, the main objective of this work is to find ways to construct initial segments which ensure maximum benefit from these infinite degrees of freedom.

## 3. Optimum utilization of the initial segment

It is shown in this section that the degrees of freedom of the initial segment can be used to couple an additional on-axis ray set. Consider the curve segments S_{1} and Ŝ_{2} of the two surfaces shown in Fig. 2(a). They are designed to couple a tilted parallel ray set onto a point. Because of the symmetry applied with respect to the optical axis of the lens, the lens segments Ŝ_{1} and S_{2} will couple the parallel ray set with opposite sign. That means that only one ray set is actually coupled using these two surface segments. This offers the opportunity to couple an additional ray set. One of these two identical off-axis ray sets is now replaced in Fig. 2(b) by an on-axis ray set. The constant optical path length condition can then be used to calculate the final lens portion on the second lens profile by tracing off-axis rays through Ŝ_{1}, which is illustrated in Fig. 2(c). Finally, the full set of rays on the completed lens is shown in Fig. 2(d) - coupling three ray sets with two lens profiles.

A simple optimization approach is used to establish the concept of convergence points and find them. These convergence points form the basis for the analytical solution given in Section 4. The implemented design algorithm starts with a central on-axis ray along the optical axis. By choosing points P_{1}, P_{2} and R_{1} the central lens thickness and the related optical path length are determined. Both normal vectors in P_{1} and P_{2} point in the direction of the optical axis due to the overall lens’ symmetry. This initial step is shown in Fig. 3(a). In the next step in Fig. 3(b), the focus position R_{2} of the off-axis ray set is chosen. A ray coming from R_{2} is refracted at P_{2}. A chosen optical path length is then used to determine the edge P_{3} on the first lens surface. At this point, all necessary information is provided - the focus positions and optical path lengths for on- and off-axis ray sets.

Figure 3(c) consists of two steps. First P_{3} and its correspondent normal vector are mirrored along the optical axis resulting in Q_{3} and its associated normal vector. An on-axis ray is refracted at Q_{3}, P_{4} is calculated using the constant optical path length condition. The correspondent normal vector is calculated to refract the ray towards the focus position R_{1}. A ray coming from R_{2} is refracted at P_{4}, P_{5} is calculated using the constant optical path length condition, as shown in Fig. 3(d). Figure 3(e) shows the next step in analogy to (c). These two steps in (c) and (d) are now alternately repeated until a stop criterion is reached.

For chosen initial values, steps (c) and (d) typically fail to close the lens profiles. The initial parameters, namely P_{1}, P_{2}, P_{3}, R_{1} and R_{2} are varied using unconstrained nonlinear optimization in MATLAB until the algorithm converges and results in the final point clouds with correspondent ray paths, which are shown in Fig. 3(f). An additional zoomed in view of one convergence point is shown to better see this behavior. The repeated design steps (c) and (d) successively construct the lens profiles from center and edge simultaneously until they asymptotically close the remaining gaps in the profiles. The optimization procedure is used to minimize the remaining gap sizes. As a consequence, in the limiting case of infinitesimal profile gaps, the design process results in two convergence points taking the axial symmetry on the first lens profile into account. The convergence points are thus characterized by the special case that on- and off-axis rays passing through the lens share identical points and normal vectors on each surface.

## 4. Analytic solution of initial segments starting from convergence points

The numerical algorithm introduced in the previous subsection provides already a comprehensive solution to the stated initial segment problem. However, multi-parametric optimization is still needed which has several disadvantages. The guess of an initial parameter set may decide whether the algorithm will converge or not, furthermore the optimization process can be fairly time consuming, to name just two.

The reverse case forms the starting point for all further considerations: determine the convergence points first and derive the full lens profiles from these initial convergence points for one on-axis and one off-axis ray set. All necessary initial values are defined as shown in Fig. 4(a).

The point coordinates (*x*_{0}, *z*_{0}) of the convergence point on the first lens profile can be freely chosen without loss of generality. The slope *m*_{0} at the convergence point represents a first variable. The intersection of the refracted on-axis ray through (*x*_{0}, *z*_{0}) and the refracted off-axis ray through the mirrored convergence point (−*x*_{0}, *z*_{0}) determines the coordinates (*x*_{1}, *z*_{1}) of the convergence point on the 2^{nd} lens profile. The slope *m*_{1} at this second convergence point represents a second variable. The intersection of the refracted on-axis ray through (*x*_{1}, *z*_{1}) and the optical axis determines the first focus position *d*. Finally, the refracted off-axis ray and the known lateral detector position determine the x-coordinate *r* of the second focus position. This formulation provides a compact representation of all initial values depending upon the two variables *m*_{0} and *m*_{1} only. In addition, it provides the permissible range of values for these two variables for which the convergence point construction is possible. For example, a given value *m*_{0} limits the possible range of values of *m*_{1} for which an on-axis ray can be refracted towards the optical axis.

Next, two analytic functions *f*(*x*) and *g*(*x*) are introduced to describe the two lens profiles. To analytically describe the optical paths of rays passing through the lens it is necessary to introduce two additional mapping functions *s*(*x*) and *u*(*x*). Figure 4(b) shows an on-axis ray passing through an arbitrary point *p⃗* = (*x*, *f*(*x*)) on the first lens profile which is then refracted in (*s*(*x*),*g*(*s*(*x*))) towards the focal point (0,*d*). The auxiliary function *s*(*x*) thus defines the mapping in x-coordinate. Similarly, function *u*(*x*) defines the mapping in x-coordinate for off-axis rays through an arbitrary point
$\widehat{\overrightarrow{p}}=(-x,f(x))$ on the first, and through (*u*(*x*),*g*(*u*(*x*))) on the second lens profile, as shown in Fig. 4(c). All optical path lengths can then be expressed in sections using vector geometry as

*n⃗*

_{0}and

*n⃗*

_{1}denote the directional vectors for on- and off-axis ray sets, respectively. The position vector

*w⃗*

_{0}denotes a fixed point on both plane wave fronts and

*n*

_{2}denotes the refractive index of the lens.

Fermat’s principle states that the optical path length between two fixed points is an extremum along a light ray. Consider a fixed point on the on-axis wave front defined by *w⃗*_{0} and *n⃗*_{0}, and the fixed point (*s*(*x*), *f*(*s*(*x*)) on the second lens profile: an on-axis ray coming from the wave front and passing through (*s*(*x*), *f*(*s*(*x*)) must be such that the combined optical path length *d*_{1} + *d*_{2} is an extremum. With point (*s*(*x*), *f*(*s*(*x*)) kept fixed, the only remaining variable to achieve an extremum for *D*_{1} = *d*_{1} + *d*_{2} is the point (*x*, *f*(*x*)) on the upper lens profile. Fermat’s principle thus implies that

*s,g*(

*s*)) is kept fixed. Similarly, an on-axis ray between fixed points (

*x*,

*f*(

*x*)) and (0,

*d*) must satisfy equation These two functional differential equations arising from Fermat’s principle describe all on-axis ray paths through the lens profiles

*f*(

*x*) and

*g*(

*x*). In analogy,

*f*(

*x*) and

*g*(

*x*) must also satisfy the functional differential equations for off-axis rays, using the same arguments as before. The lens design shown in Fig. 4 is thus fully described by four functional differential equations (Eqs. (4)–(7)) for four unknown functions

*f*(

*x*),

*g*(

*x*),

*s*(

*x*) and

*u*(

*x*). The fundamental analysis of a similar system of functional differential equations has been discussed [6–10]. The results on existence, uniqueness and analyticity were all based on fixed point theorems in functional analysis. These arguments break down when there are off-axis foci. However, the results presented in this work will provide strong evidence that the solutions for off-axis foci are analytic and smooth as well.

Suppose that (*f,g,s,u*) is an analytic and smooth solution to the functional differential equations (4)–(7), Taylor’s theorem implies that the functions must be infinitely differentiable and have a power-series representation. Thus the four functions can be given by power-series

*x*

_{0},

*z*

_{0}) and (

*x*

_{1},

*z*

_{1}), respectively. The initial conditions

*D*= 0 for

_{i}*i*= 1..4 and provide general solutions for the initial Taylor coefficients

*f*

_{0},

*f*

_{1},

*g*

_{0},

*g*

_{1},

*s*

_{0}and

*u*

_{0}depending upon variables

*m*

_{0}and

*m*

_{1}. In ascending order it is now possible to calculate (n+1)

^{th}order Taylor series coefficients in

*f*(

*x*) and

*g*(

*x*) and n

^{th}order in

*s*(

*x*) and

*u*(

*x*) by solving equations

*n*= 0 corresponds to the just solved equations for initial Taylor coefficients. There are two further cases needed to be solved:

- For
*n*= 1, the set of Eq. (11) results in nonlinear algebraic equations for Taylor series coefficients*f*_{2},*g*_{2},*s*_{1}and*u*_{1}. These equations have two general solutions, where one solution can be discarded as non-physical. The remaining unique solution can be expressed as functions of the initial, already known Taylor coefficients for*n*= 0. - For
*n*> 1, the set of Eq. (11) results in linear algebraic equations for particular Taylor series coefficients. By sorting and combining terms, the equations (11) can be transformed and expressed as a compact matrix equation$$M\left(\begin{array}{c}{f}_{n+1}\\ {g}_{n+1}\\ {s}_{n}\\ {u}_{n}\end{array}\right)=\left(\begin{array}{c}{b}_{1}^{(n)}\\ {b}_{2}^{(n)}\\ {b}_{3}^{(n)}\\ {b}_{4}^{(n)}\end{array}\right)$$for arbitrary*n*> 1. The matrix elements*M*consist of mathematical expressions which depend on Taylor series coefficients obtained for_{ij}*n*= 0,1. The vector elements of*b⃗*^{(n)}on the right hand side of Eq. (12) are mathematical expressions only dependent on previous Taylor series coefficients*f*,_{i}*g*,_{i}*s*_{i}_{−1}and*u*_{i−1}for*i*= 0..*n*. Finally, the vector elements of*b⃗*^{(n)}can be calculated for each*n*= 2,3,4,.. (in this order) from Eq. (11). For known matrix*M*and vectors*b⃗*^{(2)}..*b⃗*^{(}^{n}^{)}, the Taylor series coefficients*f*_{n+1},*g*_{n+1},*s*and_{n}*u*can then be calculated by solving the linear system of Eq. (12)._{n}

*n*= 0 and

*n*= 1 and the introduced algebraic system of linear equations (12) allow to calculate all Taylor series coefficients of

*f*(

*x*),

*g*(

*x*),

*s*(

*x*) and

*u*(

*x*) up to an arbitrary order. However, a Taylor series is a representation of an analytic function as an infinite sum of terms. Indeed, it is only possible to calculate a finite number of initial terms of the Taylor series. Such a function is called a Taylor polynomial and will be the only approximation made. Furthermore, Taylor’s remainder theorem provides quantitative estimates on convergence and the approximation error of the function by its Taylor polynomial. The radii of convergence for the expansions

*f*(

*x*) and

*g*(

*x*) are very important, as they indicate the maximum aperture that can be achieved for a given set of initial values. In the examples considered, the radius of convergence is larger than the range of 0 <

*x*<

*x*of the lenses. The algebraic steps of calculation presented in the this section are fully implemented and calculated in Wolfram

_{max}*Mathematica*.

## 5. Ray tracing results for calculated 15^{th} order Taylor polynomials

In a first step, the general solutions for Taylor polynomial coefficients for *n* = 0,1 are calculated which also determines the matrix *M* of Eq. (12). The general solution of the matrix equation *Mx⃗* =*b⃗* is then calculated to obtain the general solution vector *x⃗*. The vector elements of *b⃗*^{(n)} are calculated within a loop from *n* = 2 up to *n* = 14. It is possible to calculate even higher orders. However, the error made in this approximation is already extremely small.

All calculated mathematical expressions, sorted in the right order, are then exported as C++ compatible code and embedded in a MATLAB-compatible mex file library. Once compiled, this library returns the Taylor polynomial coefficients for the lens profiles *f*(*x*) and *g*(*x*) up to 15^{th} order, and the mapping functions *s*(*x*) and *u*(*x*) up to 14^{th} order as output arguments. The off-axis design angle *θ*, the refractive index *n*_{2} and the derivatives at the convergence points *m*_{1} and *m*_{2} are passed as input arguments.

By fixing the design angle and refractive index (here: *θ* = 12°, *n*_{2} = 1.5), the only remaining free parameters to vary are *m*_{0} and *m*_{1}. Furthermore, the intended real focus of the lens suggest to choose a negative value for *m*_{0}. The sign of *m*_{1} then determines the global lens shape - which can be either meniscus like (negative *m*_{1}) or biconvex (positive *m*_{1}). Local smooth and analytic solutions exist for any values within the permissible range of values for these two variables. However, the lens’ smoothness and symmetry additionally requires *f*′(0) = *g*′(0) = 0 at the optical axis. This boundary condition introduces an additional correlation between *m*_{0} and *m*_{1}, meaning that for a specific value *m*_{1} the correspondent value *m*_{0} can be obtained by making sure *f*′(0) = *g*′(0) = 0 is fulfilled. Even though it is not possible to deduce this correlation *m*_{1}(*m*_{0}) as a closed form solution, the numerical calculation of roots of *f*′(0) and *g*′(0) is very accurate and fast.

A ray tracing animation video (
Media 1) helps to better illustrate the correlation between *m*_{0} and *m*_{1}. For *m*_{1} ranging from −0.065 to 0.065 the correspondent *m*_{0}, satisfying the boundary condition, is directly calculated each time. Figure 5 shows two single-frame excerpts from this ray tracing animation video for *m*_{1} = −0.065 which is equivalent to a meniscus lens (a), and *m*_{1} = 0.065 which is equivalent to a biconvex lens (b). The two arrows indicate the directions of the normal vectors at the convergence points.

The ray tracing simulations for this animation video (240 frames) are done using the MATLAB-based ray tracer OPS (by courtesy of Prof. Dr. Udo Rohlfing, Hochschule Darmstadt, Germany). The ray tracer is only used for visualization purposes - the entire optics design is based on the derived general analytical solution. In addition, the two mapping functions *s*(*x*) and *u*(*x*) are used to directly calculate the semi-diameters of the lens profiles. For example, solving equation *u*(*r*_{1}) = 0 provides the semi-diameter *r*_{1} of the upper lens profile. Lens specific parameters such as the effective focal length and magnification can be directly calculated for each correspondent *m*_{0} and *m*_{1}.

Some fundamental aspects of this new analytic design method deserve to be particularly emphasized at the end of this section. The fixed parameters *θ* = 12° and *n*_{2} = 1.5 have been selected arbitrarily to provide clear evidence that it is possible to calculate analytical high order Taylor polynomials for various values *m*_{0} and *m*_{1}. However, this calculated general solution provides much more. It solves the stated initial segment problem for any given (physically meaningful) initial parameter set (*θ*,*m*_{0}, *m*_{1}, *n*_{2}). This also means that it could be used in a hybrid technique: first calculate a highly accurate initial segment and finish the overall lens design using the SMS2D design method. This particular convergence point design was chosen to demonstrate the possibility of coupling an additional on-axis ray set, meaning that three ray sets can be coupled with two surfaces. In addition, the imposed axial symmetry of the lens is not explicitly needed by this design method. Other wave fronts besides plane and spherical wave fronts could be coupled as well. A further important step will be a solution and implementation of this new method in three dimensions and for more than two surfaces. This would allow calculating free-form optics targeting applications with wide and even asymmetric field of view.

## 6. Conclusion

Within the scope of this work, it has been shown that it is possible to couple three ray sets with two refractive surfaces forming a lens of minimum thickness. The established convergence point formalism provided the basis for an analytical description of the entire lens using Fermat’s principle. The derived set of functional differential equations led to algebraic systems of equations which have been solved up to an arbitrary order of all Taylor series coefficients needed to describe the lens profiles. Exemplary ray tracing results for analytic lens profiles given as Taylor polynomials of 15^{th} order demonstrated the capabilities and versatility of this new analytic optics design method. Future work will focus on its non-rotational symmetric three-dimensional generalization and target applications where this new design method can help to further increase the overall optical system’s performance.

## Acknowledgments

Our work reported in this paper was supported in part by the Research Foundation - Flanders (FWO-Vlaanderen) that provides a PhD grant (grant number FWOTM510) and provided a grant for a research period at CeDInt, Universidad Politécnica de Madrid (grant number V424711N) for Fabian Duerr, and in part by the IAP BELSPO VI-10, the Industrial Research Funding (IOF), Methusalem, VUB-GOA, and the OZR of the Vrije Universiteit Brussel.

## References and links

**1. **M. Born and E. Wolf, *Principles of Optics* (Pergamon Press, 1986).

**2. **P. Mouroulis and J. Macdonald, *Geometrical Optics and Optical Design* (Oxford University Press, USA, 1997).

**3. **J. Chaves, *Introduction to Nonimaging Optics* (CRC, 2008). [CrossRef]

**4. **D. Grabovičkić, P. Benítez, and J.C. Miñano, “Aspheric V-groove reflector design with the SMS method in two dimensions,” Opt. Express **18**, 2515–2521 (2010). [CrossRef]

**5. **J.C. Miñano, P. Benítez, W. Lin, J. Infante, F. Muñoz, and A. Santamaría, “An application of the SMS method for imaging designs,” Opt. Express **17**, 24036–24044 (2009). [CrossRef]

**6. **A. Friedman and J.B. McLeod, “Optimal design of an optical lens,” Arch. Rational Mech. Anal. **99**, 147–164 (1987). [CrossRef]

**7. **A. Friedman and J.B. McLeod, “An optical lens for focusing two pairs of points,” Arch. Rational Mech. Anal. **101**, 57–83 (1988). [CrossRef]

**8. **J.C.W. Rogers, “Existence, uniqueness and construction of the solution of a system of ordinary functional differential equations with applications to the design of a perfectly focusing symmetric lens,” IMA J. Appl. Math. **41**, 105–134 (1988). [CrossRef]

**9. **B. van Brunt and J.R. Ockendon, “A lens focusing light at two different wavelengths,” J. Math. Anal. App. **165**, 156–179 (1992). [CrossRef]

**10. **B. van Brunt, “Mathematical possibility of certain systems in geometrical optics,” J. Opt. Soc. Am. A **11**, 2905–2914 (1994). [CrossRef]