Zheng et al.

Vol. 27, No. 8 / August 2010 / J. Opt. Soc. Am. A

1791

Analytic-domain lens design
with proximate ray tracing
Nan Zheng,1 Nathan Hagen,2,* and David J. Brady2,31
Department of Physics, Duke University, Durham, North Carolina 27708, USA
Department of Electrical and Computer Engineering, Fitzpatrick Institute for Photonics,
Duke University, Durham, North Carolina 27708, USA
3E-mail: david.brady@duke.edu*Corresponding author: nhagen@optics.arizona.edu

2Received February 17, 2010; revised June 3, 2010; accepted June 5, 2010;
posted June 7, 2010 (Doc. ID 124169); published July 14, 2010
We have developed an alternative approach to optical design which operates in the analytical domain so thatan optical designer works directly with rays as analytical functions of system parameters rather than as discretely sampled polylines. This is made possible by a generalization of the proximate ray tracing techniquewhich obtains the analytical dependence of the rays at the image surface (and ray path lengths at the exit
pupil) on each system parameter. The resulting method provides an alternative direction from which to approach system optimization and supplies information which is not typically available to the system designer. In
addition, we have further expanded the procedure to allow asymmetric systems and arbitrary order of approximation, and have illustrated the performance of the method through three lens design examples. © 2010 Optical Society of America
OCIS codes: 080.2740, 080.6755, 220.3620.

1. INTRODUCTION
The design of an optical system typically proceeds withthe use of a modern computer program for performing exact ray tracing. The user first models the optical layout
and constructs a merit function, after which the program
uses ray trace data to minimize the merit function by adjusting the system parameters. The most common optimization algorithm used is damped least-squares [1], which
is fast but has the drawback that it is limited to finding
only local minima. If the starting optical configuration is
far enough from the global minimum, and the merit function is complicated enough to possess a number of local
minima (as in almost all cases of interest), then the
chance of converging on the global solution quickly approaches zero. Global optimization techniques [2] such as
simulated annealing [3] or genetic algorithms [4] can ofcourse get around this problem, but require extensive
computational resources [5]. While it is possible to minimize the computational effort needed by providing the algorithm with a starting configuration which is as close as
possible to the global solution, thereby minimizing thesearch space, the method for doing this relies heavily on
the intuition of the designer and on tools developed inclassical aberration theory (such as aberration cancellation in symmetric systems).
We present what we believe to be a new design tool
which approaches the problem from a different direction,and thus may be helpful in situations where traditional
tools get stuck: a ray tracing engine that provides the expression for the rays, to any desired polynomial order ofapproximation, as an analytical function of the system parameters. The approach is a generalization of the proximate ray tracing method developed by Hopkins [6–8]—a1084-7529/10/081791-12/$15.00

generalization which is made possible through the powerof modern computer algebra systems [9] to manipulate
large analytical equations.
In its original form, proximate ray tracing involves first
analytically calculating the formulas for transfer and refraction at each polynomial order of approximation. A ray
trace then involves calculating the numerical values of
the ray-surface intersection point at each surface, at eachorder of approximation. The procedure was thus developed as an efficient method of numerically calculating
higher-order aberration coefficients, and it requires tracing only a small set of special rays to obtain the aberration coefficient values. The generalization we present
here extends proximate ray tracing’s analytical approachto the entire design process, without numerical substitution. Thus, rather than obtaining analytical formulas foreach individual refraction and transfer step in the trace,
we obtain a single formula for the rays at the image plane(and exit pupil). We can thus obtain the aberration coefficients in analytical form, in which the functional dependence of the lens merit function on the system parameters
is retained. The aim is to show that with the aid of modern computer algebra systems, optical design problems
that are currently performed numerically can also bedone in the analytical domain, and that this can have significant advantages for understanding the design problem.
A closely related approach for analytical ray tracing
was taken by Kondo and Takeuchi [10] through the use of
matrices and the selection of a proper vector basis formodeling the nonlinear effects present in ray tracing (upto the desired order of approximation). While similar to
the approach presented here, it lacks the conceptual sim© 2010 Optical Society of America

1792

J. Opt. Soc. Am. A / Vol. 27, No. 8 / August 2010

Zheng et al.

plicity of proximate ray tracing, and thus we feel that it isa much more cumbersome tool to work with. That is, formany optical engineers, the concept of equating terms oflike order in the Taylor expansion of a nonlinear equation
is a much more familiar process than that of constructing
a nonorthogonal vector basis in which to represent a nonlinear transformation. This conceptual simplicity becomes
important for understanding what to do with the large setof polynomial terms generated by either approach.
Kondo and Takeuchi’s matrix approach was modified by
Lakshminarayanam and Varadharajan [11], and also by
Almeida [12], and adapted for use in a computer algebra
system [13], but published work has been limited to optical modeling rather than design, in that it does not treat
the merit function or the optimization procedure. Thismisses one of the central strengths of working in the analytical domain: the tractability of polynomial equations
allows the use of more robust optimization techniques.
Other attempts at analytical ray tracing have also beenmade. Walther [14,15] developed an analytical approach
which makes use of eikonals rather than rays and is thus
less accessible for many optical engineers. Kryszczyński[16] provided some tentative steps toward analytical design based on rays, but only supplies an outline of an algorithm. Although it may at first appear mathematicallycomplex, we hope to show that the proximate ray tracing
method makes the analytical approach both accessible
and practical.
In the discussion below, we first review Hopkins’
method [6–8] and show how it can be readily generalized
to arbitrary orders of approximation and to asymmetricsystems. We then show how to construct the merit function and optimizer for designing optical systems with thisapproach and present three example designs. Finally, we
conclude with a discussion of the advantages and disadvantages of this technique.

2. TRANSFER
Hopkins [6–8] described proximate ray tracing as an iterative ray tracing technique in algebraic form. The basic
approach is
1. A polynomial series expansion for basic and intermediate variables are inserted into the exact ray trace
equations.
2. Any sines and cosines are series-expanded, and any
multiplications, divisions, and square roots are performed
as series operations.
3. Terms of a given order are collected together, andhigher-order terms are obtained using lower-order solutions via a triangular set of equations.
Ray tracing consists of two basic operations: transfer
and refraction. Using the ray path length w as the transfer parameter, the transfer equations can be written as
rs+1 = rs + wscs ,共1兲

where r = 共x , y , z兲 is the ray position vector, c = 共cx , cy , cz兲 isthe direction cosine vector, and s indicates the surface index (i.e., transfer from surface s to surface s + 1). To simplify the equations below, we will usually leave the surface index implied, except where it is needed. Here z istaken to be the optical axis, and multiplying w by the refractive index n of the medium gives the ray optical path
length. The refraction equations can be written as [[17], p.
133]
n共c ⫻ N兲 = n⬘共c⬘ ⫻ N兲,

共2兲

where N = 共Nx , Ny , Nz兲 is the normal vector of the surface,andN = ⵱f共r兲,共3兲

when f共r兲 = 0 defines the surface. Note that the optical
path length of a ray through the entire system is given by
S−1
nsws when ws is the ray path length for transfer
W = 兺s=1
from surface s to s + 1.
The first step in the proximate ray trace procedure is toexpand all of the relevant variables in the transfer and refraction equations in various orders of approximation.
Thus, each ray trace variable is expressed in the form
x = 0 + x 共1兲 + x 共2兲 + x 共3兲 + x 共4兲 + ¯ ,y = 0 + y 共1兲 + y 共2兲 + y 共3兲 + y 共4兲 + ¯ ,z = z 共0兲 + z 共1兲 + z 共2兲 + z 共3兲 + z 共4兲 + ¯ ,cx = 0 + cx共1兲 + cx共2兲 + cx共3兲 + cx共4兲 + ¯ ,cy = 0 + cy共1兲 + cy共2兲 + cy共3兲 + cy共4兲 + ¯ ,cz = 1 + cz共1兲 + cz共2兲 + cz共3兲 + cz共4兲 + ¯ ,w = w 共0兲 + w 共1兲 + w 共2兲 + w 共3兲 + w 共4兲 + ¯ .

共4兲The superscripts in parentheses indicate the order of approximation so that the first nonzero term on the righthand side of each of these equations represents a paraxial
variable, and succeeding terms represent the nonlinear
dependence on the paraxial variables. The term z共0兲 represents the axial transfer distance from the surface vertex
to the previous surface (and is thus a negative quantity
for rays propagating from left to right).
The primary variables used to define the rays, the entrance pupil coordinates 共xep , yep兲 and the field angles
共␪x , ␪y兲, are treated as paraxial variables and thus do nothave an order-expansion. The final expressions for therays will give the image coordinates in terms of these primary variables and of the parameters used to define eachsurface. In aberration theory it is more common to workwith normalized field angles H, defined either as
共Hx , Hy兲 = 共1 / ␪max兲共␪x , ␪y兲,where␪max ⬅ 关共max兵␪x其兲2+ 共max兵␪y其兲2兴1/2, or as 共Hx , Hy兲 = 共␪x / max兵␪x其 , ␪y / max兵␪y其兲.
In the presentation below, we will continue to use ␪ rather
than H to represent field angles.
The elements of the direction cosine vector for the incident ray are given by cx = sin ␪x, cy = sin ␪y, and cz= ± 冑1 − c2x − c2y for a ray propagating in the ±z direction.
Thus
共cx共1兲 + cx共3兲 + ¯兲 = ␪x −13!

␪3x + ¯ ,Zheng et al.

Vol. 27, No. 8 / August 2010 / J. Opt. Soc. Am. A

共cy共1兲 + cy共3兲 + ¯兲 = ␪y −13!

␪3y + ¯ ,

冋1
1
1
共␪2x + ␪2y 兲2 + 共␪4x
共cz共0兲 + cz共2兲 + cz共4兲 + ¯兲 = 1 − 共␪2x + ␪2y 兲 +
2
16
6

册

+ ␪4y 兲 + ¯ .
In Hopkins’ presentation [6–8] of the proximate ray
trace equations, the variables x, y, cx, and cy have only
odd-order terms, and the variables z, cz, and w have only
even-order terms due to his assumption of rotational symmetry about the optical axis. In order to design more general optical systems, we drop these symmetry assumptions here and use the most general form of the equations.
Substituting the order-expanded variables for r, c, andw into the transfer equations (1) results in
0 = z0共0兲 + w共0兲 ,1793

the first-order equation allows one to solve for w共1兲, andwe can likewise continue to substitute lower-order results
to obtain solutions for the higher-order equations. The resulting sequence of values for w, i.e., 共w共0兲 , w共0兲 + w共1兲 , w共0兲
+ w共1兲 + w共2兲 , . . .兲, provides an estimate of the real ray path
length to increasing order of approximation. Once all desired orders of w have been solved for, one can then substitute into the equations for x and y. As with w, the se共0兲
共1兲 共0兲
共1兲
共2兲
quence 共r共0兲
provides ans , rs + rs , rs + rs + rs , . . .兲estimate of the ray-surface intersection location to increasing order of approximation (see Fig. 1).
One further step is necessary before we can use thisprocedure to solve this set of equations. Since the various
orders of the surface sag zs are not yet known, we cannot
yet solve directly for w. First we need to express zs interms of known quantities, and for this we need to perform the order-expansion of the surface equation.

3. SURFACE EQUATIONThe order-expansion of the surface equation is obtainedby taking its Taylor expansion, shown here for a surface
in the form z = z共x , y兲

r共s1兲 = r0共1兲 + c共1兲w共0兲 + c共0兲w共1兲 ,冋 册 冋 册冋 册 冋 册冋 册

r共s2兲 = r0共2兲 + c共2兲w共0兲 + c共1兲w共1兲 + c共0兲w共2兲 ,⳵ zs

zs共xs,ys兲 = zs共0,0兲 + xs

r共s3兲 = r0共3兲 + c共3兲w共0兲 + c共2兲w共1兲 + c共1兲w共2兲 + c共0兲w共3兲 ,

⳵ xs

1
⳵ zs
+ x2s
2
⳵ x2s

in which we label the surface before transfer as s = 0 andthat following transfer as simply s. Each line contains
only terms of equal order, such that the number of equations in the resulting system is equal to the desired order共0兲
of approximation. Using the fact that cz共0兲 = 1 and c共0兲
x = cy= 0, the transfer equations split into

共0,0兲

1
⳵ 2z s
+ y2s
2
⳵ y2s

共0,0兲

w共0兲 = − z0共0兲 ,w共1兲 = z共s1兲 − z0共1兲 − cz共1兲w共0兲 ,

⳵ ys

共0,0兲

共0,0兲

⳵ zs

2

]

⳵ zs

+ ys

2

+ x sy s

⳵ xs ⳵ ys

共0,0兲

+ ¯,

共6兲

in which the 共xs , ys兲 = 共0 , 0兲 subscript on each square
bracket indicates that the partial derivatives are evaluated at the axial surface point. From Eq. (6), we next substitute in the order-expansion forms of 共xs , ys , zs兲 andequate terms of equal order, giving
z共s0兲 = zs共0,0兲,

w共2兲 = z共s2兲 − z0共2兲 − cz共2兲w共0兲 − cz共1兲w共1兲 ,]z共s1兲 = x共s1兲

x共s1兲 = x0共1兲 + cx共1兲w共0兲 ,x共s2兲 = x0共2兲 + cx共2兲w共0兲 + cx共1兲w共1兲 ,

冋 册
⳵ zs

⳵ xs

+ y共s1兲
共0,0兲

( y ,z

]

⳵ zs

⳵ ys

,
共0,0兲

(y'(1),z'(0))( y (1),z(0))
(3)

冋 册

(y'(3),z'(1))

(1))

t

y共s1兲 = y0共1兲 + cy共1兲w共0兲 ,y共s2兲 = y0共2兲 + cy共2兲w共0兲 + cy共1兲w共1兲 ,]共5兲

The zeroth-order equation gives w共0兲 = −z0共0兲, which is simply the axial distance from the previous surface vertex tothe current surface vertex. Substituting this result into

Fig. 1. (Color online) The transfer operation involves taking
low-order polynomial approximations of surfaces and modifying
the surface intersection coordinates to increasing accuracy as
higher orders are traced. Shown here are only the first- andthird-order approximations of a spherical surface. The ray path
length w from the left to right surfaces depends on the order ofapproximation, as indicated in Eqs. (5).

1794

J. Opt. Soc. Am. A / Vol. 27, No. 8 / August 2010

z共s2兲 = x共s2兲冋 册 冋 册冋 册冋 册 冋 册
⳵ zs

⳵ xs

⳵ zs

+ y共s2兲

⳵ ys

共0,0兲

1
⳵ zs
+ 关y共s1兲兴2
2
⳵ y2s

1
⳵ 2z s+ 关x共s1兲兴2
2
⳵ x2s
共0,0兲

⳵ zs

2

Zheng et al.

the normalization condition for the direction cosine vector, 兩c⬘兩 = 1, we can solve for the refracted direction cosines
c ⬘,

共0,0兲

2

+ x共s1兲y共s1兲

共0,0兲

⳵ xs ⳵ ys

,

cx⬘ =

共0,0兲

]

cy⬘ =

In the case of a spherical surface with radius of curvature
R whose center of curvature lies on the optical axis, thesurface equation can be written as z = −R + 冑R2 − r2s (for rs= 冑x2s + y2s ). The order-expansion of this surface producesz共s0兲 = 0,r共s1兲r共s3兲
R

+

cz⬘ =

Axy
Bz ± 冑D
Az

共8兲

,

共9兲

,

共10兲

Az = n22共N2x + N2y + Nz2兲,,
Bx = n1n2Nz关N2y cx − NxNycy + Nz共Nzcx − Nxcz兲兴,z共s3兲 = 0,z共s4兲 =

By ± Ny冑D,

Axy = n22Nz共N2x + N2y + Nz2兲,关r共s1兲兴2
2R

Axy

where

z共s1兲 = 0,z共s2兲 =

Bx ± Nx冑D

By = n1n2Nz关N2x cy − NxNycx + Nz共Nzcy − Nycz兲兴,关r共s1兲兴4
8R

Bz = n1n2关− Nz共Nxcx + Nycy兲 + 共N2x + N2y 兲cz兴,,3

D = n22Nz2关n22共N2x + N2y + Nz2兲 − n12共Nz2共c2x + c2y 兲 − 2NxNzcxcz
]
共1兲 2
共1兲 2
2
where 关r共1兲
s 兴 ⬅ 关xs 兴 + 关ys 兴 . The even-order components
are equal to zero here due to rotational symmetry. Now
we can see that higher-order forms of zs can be written interms of lower-order forms of xs and ys. Thus, if we returnto Eqs. (5) and substitute in for the order-expanded formof zs given here, we find that the higher-order equations
can all be expressed in terms of lower-order quantities, allowing the full system of equations to be solved.
While we have shown the order-expansion of a spherical surface, one may also define a surface by an arbitrary
polynomial in xs and ys, such that z = z共x , y兲 is given byzs = ␣1xs + ␣2ys + ␣3x2s + ␣4xsys + ␣5y2s + ␣6x3s + ␣7x2s ys+ ␣8xsy2s + ␣9y3s + ¯ .

共7兲

In this case the surface order-expansion gives− 2Nycy共Nxcx + Nzcz兲 + N2y 共c2x + cz2兲 + N2x 共c2y + cz2兲兲兴,共11兲
and n1, n2 are the refractive indices of the media before
and after refraction. In the equations for c⬘ [Eqs.
(8)–(10)], choosing for the solution the positive sign infront of the square root selects a ray propagating in the +z
direction.
The square root, multiplication, and division in Eqs.
(8)–(10) are each nonlinear procedures and so we must
perform each operation in the context of power series tothe appropriate order. The order-expansion of the squareroot can be done by searching for an order-expanded variable ␣, whose square is equal to D, i.e.,
共␣共0兲 + ␣共1兲 + ¯兲共␣共0兲 + ␣共1兲 + ¯兲 = 共D共0兲 + D共1兲 + ¯兲.
This involves solving a triangular set of equations,z共s0兲 = 0,

␣ 共0兲␣ 共0兲 = D 共0兲 ,z共s1兲 = ␣1x共s1兲 + ␣2y共s1兲 ,␣ 共0兲␣ 共1兲 + ␣ 共1兲␣ 共0兲 = D 共1兲 ,z共s2兲 = ␣1x共s2兲 + ␣2y共s2兲 + 2␣3关x共s1兲兴2 + ␣4x共s1兲y共s1兲 + 2␣5关y共s1兲兴2 ,␣ 共0兲␣ 共2兲 + ␣ 共1兲␣ 共1兲 + ␣ 共2兲␣ 共0兲 = D 共2兲 ,]]The next step in the ray trace procedure is solving for therefracted ray direction. By combining the refraction equations (2) and the equation for the surface normal (3) with

Substituting the order-expanded variables into definition(11) of D and sorting terms by order, we can obtain theexpressions for D共0兲, D共1兲, etc. Inserting these into Eqs.
(12), we can solve the zeroth-order equation to give ␣共0兲.
Following the back-substitution procedure, we then use
each lower-order solution to solve each higher-order equa-

4. REFRACTION

共12兲

Zheng et al.

Vol. 27, No. 8 / August 2010 / J. Opt. Soc. Am. A

1795

tion and eventually obtain all of the unknown terms in ␣
up to the desired order.
The next step is to perform the division step in Eqs.
(8)–(10). For a division such as c = b / a in which all threevariables are taken to have series form (i.e., c = c共0兲 + c共1兲+ c共2兲 + ¯, etc.), we can solve this operation by first multiplying both sides by a,so that the division, square root, and square operations
must be done in order-expansion form. In the case of thepolynomial surface example (7), however, we readily obtain the result explicitly,

共c共0兲 + c共1兲 + ¯兲共a共0兲 + a共1兲 + ¯兲 = 共b共0兲 + b共1兲 + ¯兲,Nx共0兲 = − ␣1 ,and once again solving the resulting triangular set of
equations for the terms of c,

Ny共0兲 = − ␣2 ,c 共0兲a 共0兲 = b 共0兲 ,Nx共1兲 = − 2␣3x共s1兲 − ␣4y共s1兲 ,

⳵
⳵y

−y
z共x,y兲 =冑R − x 2 − y 22,

Ny共1兲 = − ␣4x共s1兲 − 2␣5y共s1兲 ,c 共0兲a 共1兲 + c 共1兲a 共0兲 = b 共1兲 ,Nx共2兲 = − 2␣3x共s2兲 − ␣4y共s2兲 − 3␣6关x共s1兲兴2 − 2␣7x共s1兲y共s1兲 − ␣8关y共s1兲兴2 ,c 共0兲a 共2兲 + c 共1兲a 共1兲 + c 共2兲a 共0兲 = b 共2兲 ,Ny共2兲 = − ␣4x共s2兲 − 2␣5y共s2兲 − ␣7关x共s1兲兴2 − 2␣8x共s1兲y共s1兲 − 3␣9关y共s1兲兴2 ,]]via back-substitution.
Finally, we also need to obtain the equation for the surface normal vector N = ⵱关z共x , y兲 − z兴 so that
Nx =

Ny =

⳵
⳵ xs
⳵
⳵ ys

zs共xs,ys兲,

zs共xs,ys兲,

共13兲Nz = − 1.

As with all other nonprimary variables in the system, wetake an order-expansion of the normal vector components
in terms of primary variables. The order-expansion for N
takes the form

Performing this sequence of operations once for each ofEqs. (8)–(10) gives the solution for the refracted ray direction cosine vector c⬘. The resulting sequence of values forc⬘, i.e., 共c⬘共0兲 , c⬘共0兲 + c⬘共1兲 , c⬘共0兲 + c⬘共1兲 + c⬘共2兲 , . . .兲, provides anestimate of the real refracted ray angle to increasing order of approximation. Note that for an nth-order ray
trace, the surface must be expanded to order n + 1 prior to
taking its derivative in order to obtain an nth-order form
for the surface normal.
While following this ray trace procedure manually is tedious and error-prone, it can be made fast and robust
through the use of modern computer algebra systems. In
fact, most such systems provide enough functionality thatthe entire transfer and refraction operations can each beperformed in a couple lines of code, and is typically executed within seconds for systems of modest complexity.
(See Appendix A for comments on how to structure thecode to help make this possible.)

Nx = Nx共0兲 + Nx共1兲 + Nx共2兲 + Nx共3兲 + Nx共4兲 + ¯ ,5. MERIT FUNCTION

Ny = Ny共0兲 + Ny共1兲 + Ny共2兲 + Ny共3兲 + Ny共4兲 + ¯ ,Nz = Nz共0兲 + Nz共1兲 + Nz共2兲 + Nz共3兲 + Nz共4兲 + ¯ ,共14兲which we can use to replace the terms on the left handside of each equation in Eqs. (13). On the right hand sideof each equation, we can use the computer algebra system
to obtain the derivative of zs and substitute into the result
the order-expansion forms of xs and ys. In general, this requires a great deal of analytical work to perform eachmathematical operation in order-expansion form, but
computer algebra systems can work through these steps
without difficulty. For example, for a spherical surface,
the partial derivatives needed for Eqs. (13) are

⳵
⳵x

−x
z共x,y兲 =冑R 2 − x 2 − y 2,

The final result of a proximate ray trace calculation is toobtain x = x共1兲 + x共2兲 + ¯ and y = y共1兲 + y共2兲 + ¯, the position of
the ray at the image plane to each order of approximation.
Each term is itself a function of the ray coordinate at theentrance pupil 共xep , yep兲 and the incident ray angle 共␪x , ␪y兲
so that we have a polynomial expression in these four
variables in addition to all of the parameters used to define the various surfaces and their spacings. In order todesign a system, we need to construct a merit function,
for which the mean square spot size is a common choice[18]. Denoting the merit function by M, we can write
M=

冕冕冕冕

关x共 兲 − xG共 兲兴2 + 关y共 兲 − yG共 兲兴2dxepdyepd␪xd␪y ,共15兲

where 共 兲 represents the variable and parameter dependence of the ray coordinates, i.e., 共xep , yep , ␪x , ␪y , . . .兲, in

1796

J. Opt. Soc. Am. A / Vol. 27, No. 8 / August 2010

which the ellipsis indicates the surface parameters of thesystem. For a system of spherical surfaces in a rotationally symmetric design, the surface parameters take theform 共R0 , t0 , n0 , R1 , t1 , n1. . .兲—the radius of curvature,
spacing, and refractive index for the successive surfaces.
In the merit function integral, 共xG , yG兲 is the Gaussian image point, for which we can simply substitute the firstorder solution of the ray location at the image plane,
共x共1兲 , y共1兲兲.
As in the ray trace procedure, if we wish to work withan analytical merit function, we can insert the orderexpanded variables for x共 兲 and y共 兲 in the above integral,collect terms of like order, and perform the integral oneach order independently. The resulting expression can bequite long for optical systems of even moderate complexity, and so the use of a computer algebra program is essential here. If the optical system possesses rotational
symmetry, we can simplify the merit function expression
to have the formM=冕冕冕关x共 兲 − xG共 兲兴2 + 关y共 兲 − yG共 兲兴2d␾␳d␳d␪ ,in which the primary variables are no longer
共xep , yep , ␪x , ␪y , . . .兲 but rather 共␳ , ␾ , ␪ , . . .兲. That is, the field
angle can be expressed as a scalar, and the pupil location
is now expressed in cylindrical coordinates, i.e., 共xep , yep兲
= 共␳ cos ␾ , ␳ sin ␾兲. For the majority of imaging systems,
the most appropriate choice of integration range is a rectangular field and a circular pupil. As long as the integration range allows us to obtain analytical functions forpolynomial integrands, then it remains possible to obtainan analytical function for the merit function as well. Notethat the squaring of the terms in the integrand results ina merit function polynomial of order 2p + 4 after integration, where p is the order of approximation in the ray
trace, and the additional four orders arise from the four
integrals of Eq. (15).
In addition, if we wish to use the spot centroid rather
than the Gaussian image point as our reference for the
merit function, a choice which amounts to ignoring the effects of distortion on the image, then we can replace
共xG , yG兲 with the appropriate centroid 共x̄ , ȳ兲 given byx̄共 兲 =

冕冕冕冕x共 兲dxepdyepd␪xd␪y ,ȳ共 兲 =

冕冕冕冕y共 兲dxepdyepd␪xd␪y .

After constructing the merit function, the final step indesigning an optical system is the implementation of an
optimization algorithm to determine the system parameters which minimize M. Here we run into many of thesame problems encountered by optimization in the existing design software: while local techniques are compact
and fast, they typically cannot reach the global solution;
while global techniques are capable of finding the optimal
solution, they require unrealistic computational resourcesin order to do so. (Reference [19] provides a useful surveyof modern algorithms for solving polynomial equations.)
For rotationally symmetric systems of modest complexity

Zheng et al.

or asymmetric systems of low complexity, existing algorithms are capable of locating global minima. Beyond
these, one must compromise between the computational
resources available and the restriction to local domains.
In the examples shown in Section 7 below, for low complexity systems, we use fast global techniques such as
Mathematica’s [20] NMinimize function, whereas formore complex problems we resort to simulated annealing.
If the designer wishes, it is also possible to constrain
the optimization using implicit functions of the systemparameters. For example, if we wish to restrict the lensdiameters to be within some allowed range, then—given
the functional form of the ray at the appropriate
surface—we can obtain equations of constraint. For example, for a rotationally symmetric spherical lens, thesurface equation gives
y2
z共y兲 =

2R

y4
+

8R3

+ ¯,

for ray height y and radius of curvature R. Constraining y
to be less than some value y0 allows us to solve for an
equation of constraint on R. This can be used by the optimization routine to look for solutions lying only within
the valid design space.
The optimization can run into trouble due to the sheer
size of the analytical formulas produced by the ray trace,
especially for systems with more than a few surfaces andwith surfaces having many modeling parameters (such as
high-order aspheres). When this happens, one thing thatcan be done is to fix some of the system parameters andoptimize over the remaining ones. For example, if we fix
the thickness of a lens, then we can give it the numericalvalue during ray tracing so that it need not be tracked
analytically. This can greatly simplify the resulting expressions and make ray tracing and optimization much
faster.

6. EXAMPLE RAY TRACE
Since the discussion up to this point has given the general
expressions for the proximate ray tracing technique, we
illustrate the approach with a simple example, using thelens shown in Fig. 2 (see also Table 1). This lens model
has been chosen such that the example ray trace can bepresented easily on a printed page: an f / 2.2 singlet cylin-

Fig. 2. (Color online) 2D lens model used for the example. Oncethe analytical model is completed, we substitute the following parameter values to give the lens shown: t1 = 5 mm, t2= 76.6667 mm, and the refractive index is n = 1.5. The first lenssurface is convex, and the second surface is planar. The resulting
f / 2.2 design has an entrance pupil diameter of 32 mm and a focal
length of 80 mm.

Zheng et al.

Vol. 27, No. 8 / August 2010 / J. Opt. Soc. Am. A

Table 1. Prescription for the Lens Shown in Fig. 2

Pupil
Lens front
Lens backImages

Radius of CurvatureThickness

Index

0
1
2
3

⬁
R
⬁
⬁

0
t1
t2

n

2R2

2
2
关t1yep
共n − 1兲2 + ␪2y R2t1n2 − yepR + 2␪yyepRt1n共ny2共1兲 =y2
2R

1

− 1兲兴,

drical lens with the (planar) aperture stop placed to coincide with the front surface of the lens. Limiting the ray
trace to two dimensions reduces the definition of a ray toits 共y , z兲-coordinates and its angle ␪y. In the ray trace
given below, all subscripts refer to the surface number s,
and we limit the approximation to third order for spacereasons. To match common practice in ray tracing, thesurface equations are expressed in terms of the surface
vertex point, with the axial transfer distance inserted into
w共0兲 in the transfer equations.
The ray trace starts at the pupil plane 共s = 0兲, with theprimary variables yep and ␪y defining all incident rays.
The first surface is defined by the equationz=

w1共2兲 =

y2共3兲 =

1
6R3

8R3

t1
R

yep共n − 1兲 + yep + ␪yt1n,

共3␪2y yepR2t1n共n − 1兲共1 + 3n兲 + ␪3y R3t1n共3n2 − 1兲3
2
共n − 1兲关− R + t1 + t1n共n − 1兲兴 − 3␪yyep
R共n − 1兲关R
+ 3yep+ t1共3n2 + n − 1兲兴兲,
z2共0兲 = 0,z2共2兲 = 0,
and again refracting,

y4
+

1797

.

关cy共1兲兴2 =1
6nR3

关6yepR2共n − 1兲 − 6␪yR3n兴,

Applying the transfer equations gives
w0共0兲 = 0,

w0共2兲 =

y1共1兲 = yep,

y1共3兲 =

z1共0兲 = 0,z1共2兲 =

2
yep

2R

,

2
␪y
yep

2R
2
yep

2R

关cy共3兲兴2 =关cy共3兲兴1 =

1
6R3

冋

1
R

关cz共0兲兴2 = 1,

关cz共2兲兴2 =,册

共n − 1兲yep + ␪yn ,2共n2 − n − 1兲兴,− 3␪yyep关cz共0兲兴1 = 1,
关cz共2兲兴1 = −2R23关␪3y R3n − 3␪2y yepR2n共n − 1兲 − 3yepn共n − 1兲,3关␪3y R3n − 3␪2y yepR2n共n − 1兲 − 3yepn共n − 1兲

1

6R3n

2
R共2n2 + n − 1兲兴,+ 3␪yyepand the refraction equations for the refracted direction cosines are关cy共1兲兴1 =

−1

−1
2R2n2

共关yep共n − 1兲 + ␪yRn兴2兲.

Finally, we transfer to the image surface. Here the only
quantity of interest is the ray coordinate y, so we omit theexpressions for w and z,
y3共1兲 =

y3共3兲 =

1Rn

关yep共− t2 + nR − nt1 + nt2 + t1n2兲 + ␪yR共nt2 + n2t1兲兴,

1
6R3n3

3
兵yep
关3t2共n4 − 3n2 + 3n − 1兲 + 3n3共n − 1兲共− R + t1

2
+ t1n共n − 1兲兲兴 + ␪yyep
关3Rt2n共n − 1兲共2n2 + 4n − 3兲

− 3Rn3共n − 1兲共R + t1共− 3n2 + n − 1兲兲兴+ ␪2y yep关3R2t2n2共n − 1兲共n + 3兲 + 3R2n4t1共n − 1兲共3n2
关yep
共1 − 2n + n2兲 + yep2n共␪yRn − ␪yR兲+ ␪2y R2n2兴.
Transferring to the planar back surface of the lens,
w1共0兲 = t1 ,+ 1兲兴 + ␪3y 关2R3t2n3 + R3t1n4共3n2 − 1兲兴其.
This is the analytical expression giving the location of allrays at the image plane in terms of the variables defining
the incident ray 共yep , ␪y兲, and of the system parameters
共R , t1 , t2 , n兲, to third-order approximation. In order to use
this model to design a lens, we construct a merit function,

1798

M=

J. Opt. Soc. Am. A / Vol. 27, No. 8 / August 2010冕冕关y3共 兲 − y3共1兲共 兲兴2dyepd␪y =冕冕关y3共3兲共 兲兴2dyepd␪y .

At this point, we have not defined any of the system parameters, pupil size, or field of view. Defining the latter
two allows us to perform the above integrals, and a typical optical design procedure would involve fixing the values of t1, t2, and n—or constraining their ranges—prior tosearching for the optimum. (If these three parameters areallowed to vary freely, one can always obtain a zeroaberration solution by letting n → ⬁ or t2 → ⬁.) For a ⫾10°
field of view and a 32 mm pupil diameter, and fixing thesystem parameters at t1 = 5 mm, t2 = 76.6667 mm, and n= 1.5, we find the optimal radius of curvature to be R
= 41.75 mm.
Even for such a simplified case, we can see that the expressions are lengthy so that performing design in theanalytical domain requires interacting with these functions via computer algebra systems. They do, however,provide a much more informative description of the rayintercepts’ nonlinear dependence on each given systemparameter. This can create a problem of information overload, in contrast to the information uncertainty producedby sampled exact ray tracing.
Note that one needs not expand all surfaces to the samepolynomial order of approximation. The only restriction
here is that the order of approximation in the ray traceneeds to be at least as much as the highest order used forthe surfaces within the system.

7. DESIGN EXAMPLES
To illustrate the performance and flexibility of our design
approach, we show the design of three example systems
and compare to results derived with conventional opticaldesign software. While comparisons between conventional and new techniques developed during research arealmost invariably unfair, in that more effort is put into
optimizing the latter than the former for the specific cases
at hand, what we wish to show is that the approach differs qualitatively from conventional methods and is capable of providing equivalent answers in a wide variety ofdesigns so that it shows promise as a new design tool.
The first example is a variant of the example ray trace
performed in Section 6, but this time using the lens offaxis and allowing the two surfaces to be freeform x-y polynomials rather than spherical (see also Fig. 3). The sys-

Fig. 3. (Color online) The lens used in the first design example:
a freeform x-y polynomial singlet lens with an off-axis field. The
system parameters are similar to those of Fig. 2: t1 = 5 mm, t2= 76.6667 mm, and the refractive index is n = 1.5, while the surface parameters are the 14 polynomial coefficients (seven foreach surface) out to fourth order. The field of view for this designis 5 ° ⱕ ␪y ⱕ 7° and −1 ° ⱕ ␪x ⱕ +1°.

Zheng et al.

tem parameters are t1 = 5 mm, t2 = 76.6667 mm, therefractive index n = 1.5, and the rectangular field of viewis 5 ° ⱕ ␪y ⱕ 7° and −1 ° ⱕ ␪x ⱕ +1°. Thus, the equations forthe front and back surfaces of the lens are defined as
zfront = ␣1x2 + ␣2y2 + ␣3y3 + ␣4x2y + ␣5x4 + ␣6y4 + ␣7x2y2 ,zback = ␤1x2 + ␤2y2 + ␤3y3 + ␤4x2y + ␤5x4 + ␤6y4 + ␤7x2y2 ,with the ␣s and ␤s coefficients left to be determined by the
optimizer. Since the system is symmetric about the x = 0
plane, all odd-order terms in x are necessarily zero.
In order to compare the results with conventional design approaches, we designed this lens using both theproximate ray trace technique and Zemax [21], with theresulting coefficients shown in Table 2. The figure of meritfrom each design can be calculated in either the domainused by proximate ray tracing (i.e., an approximate analytical model) or the domain used by Zemax (i.e., an exact
sampled model), with the following results:
Merit Function

Design Method

DomainProximate

Zemax

Analytic Approx.
Sampled Exact

0.001180.07483

0.00205
0.05635

Note that both merit functions require approximations—
truncated order in the case of proximate ray tracing andsampling in the case of Zemax’s default method.
The second design example shows a setup appropriate
to the design of a lenslet used in a multiscale lens [22]. A
multiscale lens involves the use of a standard objective
lens combined with a back-end lenslet array used to perform remapping and aberration correction on the image
prior to detection. The modeling of these systems can bequite complex since the lenslet elements are designed tohave freeform surfaces, and pupil vignetting is both large
and varies rapidly with field angle. The number of surfaces present in the system is small, however, making the
problem tractable for analytical ray tracing.
Table 2. Surface Parameters Obtained for theFirst Design Example (an Off-Axis Singlet)a
Proximate

␣1
␣2
␣3
␣4
␣5
␣6
␣7
␤1
␤2
␤3
␤4
␤5
␤6
␤7

−2

1.188⫻ 10
1.137⫻ 10−2
3.620⫻ 10−6
1.116⫻ 10−5
−3.863⫻ 10−7
−1.573⫻ 10−6
1.839⫻ 10−6
−5.649⫻ 10−4
−9.805⫻ 10−4
2.782⫻ 10−6
1.150⫻ 10−5
−1.014⫻ 10−6
−2.030⫻ 10−6
1.034⫻ 10−6

Zemax
1.136⫻ 10−2
1.114⫻ 10−2
−2.893⫻ 10−5
−1.790⫻ 10−5
2.478⫻ 10−7
−1.097⫻ 10−6
−8.058⫻ 10−7
−1.124⫻ 10−3
−1.234⫻ 10−3
−2.856⫻ 10−5
−1.675⫻ 10−5
−1.410⫻ 10−7
−1.400⫻ 10−6
−1.467⫻ 10−6aThe proximate ray trace is performed out to sixth order and optimized withsimulated annealing; the Zemax results use damped least-squares optimization.

Zheng et al.

Vol. 27, No. 8 / August 2010 / J. Opt. Soc. Am. A

This example illustrates the design of a single lenslet
within the array. The goal of the design is to re-image asection of the initial image plane and also to perform aberration correction so that the re-imaged field contains
less blur than the original field. In the case here, the field
angles re-imaged by the lenslet are 5.75° ⬍ ␪ ⬍ 8.25°, theentrance pupil diameter is 8 mm, and the lenslet itself is
4 mm in diameter. For the lenslet shown here, the freeform rear surface is modeled using an x-y polynomial
while the front surface is spherical so that the rear surface has the equationz = ␣ 1x 2 + ␣ 2y 2 + ␣ 3y 3 + ␣ 4x 2y + ␣ 5x 4 + ␣ 6y 4 + ␣ 7x 2y 2 .
The optical layout and prescription for this setup are
given in Fig. 4 and Table 3.
Once again performing the proximate ray trace to sixthorder and constructing the merit function, the simplicityof this problem allows us to use the general-purpose
NMinimize function in Mathematica. Likewise using Zemax’s default optimization tool (damped least-squares),we obtain the following results:
Merit Function

1799

Table 3. Multiscale Lens Design Example Layout
and PrescriptionsRadius of Curvature(mm)

Thickness
(mm)

Pupil
0
⬁
Obj. lens front 1
79.60
Obj. lens back2
−55.80First image
3⬁
Coord. break: rotate about x axis by 4.54°
Lenslet front
4
7.88
Lenslet back5⬁
Final image
6
⬁

20
6
72.02
62.44
1
9

Index

1.5168

1.8

gion is thus 2 mm long, and we wish to maximize theamount of incident light reaching this region. If we confine ourselves to analyzing the system in two dimensions,then the surface equations for this lens are
zfront = ␣1y2 + ␣2y4 + ␣3y6 ,zrear = ␤1y2 + ␤2y4 + ␤3y6 .

Design Method

DomainProximate

Zemax

Analytic Approx.
Sampled Exact

0.0003428
0.127130.002169
0.01467

A naive attempt at constructing a merit function forthis problem would be something likeM=

As before, we find that each method is optimal in its own
domain, although the two designs are quite similar inshape (see Table 4).
The third design example is a rotationally symmetricnonimaging concentrator, where each of the lens surfaces
is even aspheric. The goal here is different from that of
the two previous examples in that we are no longer concerned with imaging quality per se. Rather, we want tomaximize the amount of light we can concentrate onto agiven region at the “image plane”—what we call the concentration region. Thus, for a given range of incident ray
angles, we attempt to maximize the number of rays incident on the concentration region.
Due to symmetry, we can produce an approximate solution by confining the field angles and pupil coordinates
to the meridional plane. Thus, we consider a field of viewof −20° ⱕ ␪y ⱕ +20°; the entrance pupil diameter is 2.93mm, the system track length is 3.14 mm, and the image
size is −1 mmⱕ y ⱕ +1 mm. (This design example is modeled after [[23], pp. 189–192].) The light concentration re-Fig. 4. (Color online) The multiscale lens design example layout
and prescription. The objective is fixed, while we attempt to design the lenslet to perform aberration correction on the nominal
image (shown by the curved surface between the lenses) and reimage onto a tilted detector array (shown at the far right). The
square pupil is 8 mm⫻ 8 mm in size, and the objective lens focal
length is 64.45 mm. The lenslet re-images a 2.5° ⫻ 2.5° square
field of view.

冕

Dep/2

dyep

−Dep/2冕20°d ␪ yy 2 ,−20°

where y gives the position of the ray at the image plane
and Dep is the diameter of the entrance pupil. While thisfunction penalizes rays which stray too far from the axis,
what we really want is a penalty which is zero or verysmall for rays falling onto the concentration region, butvery large for rays falling outside it. The implementationof this approach is tricky, however, as it requires ᐉ1 minimization rather than the much more widely used ᐉ2 minimization techniques, and also it requires that we work directly with image coordinates as primary variables,
rather than the object coordinates we have been using up
to now. This is a topic we hope to treat at length in a future publication.
An alternative optimization approach takes advantageof the edge ray principle [[23], p. 183]. The phase space forrays propagating through the system (Figs. 5 and 6) arebounded by the square abcd. We can choose to map those
Table 4. Surface Parameters Obtained for theSecond Design Example (Multiscale Lenslet)a

␣1
␣2
␣3
␣4
␣5
␣6
␣7

ProximateZemax−0.070 26−0.068 24
0.000 08
0.000 020.002 35
0.002 240.004 02

−0.066 25−0.064 25
0.000 190.000 170.001 380.001 310.002 38a
The proximate ray trace is performed out to sixth order and optimized with
Mathematica’s NMinimize function; the Zemax results use damped least-squares optimization.

1800

J. Opt. Soc. Am. A / Vol. 27, No. 8 / August 2010

Fig. 5.

Phase space of the incident rays.

rays represented by line aa⬘ to the edge of the concentration region—a point we designate as ym—and map therays represented by line a⬘b to the edge of the back surface of the lens—a point we designate as 共yk , zk兲. By symmetry, this likewise forces all the rays represented by line
cc⬘ to focus near −ym, and all the rays represented by linec⬘d to focus near 共−yk , zk兲. The corresponding merit function can be given the formM=冋冕关共y3 − yk兲2 + 共z3 − zk兲2兴dyep冕共Dep/2兲−ya⬘−Dep/2共y − ym兲2dyep

册

,
␪=20°

where ya⬘ is the yep value of ray a⬘. (The value of ya⬘ isdetermined during the optimization step.) The coordinate
共y3 , z3兲 is the ray position on the back surface of the lens—surface 3 in the system. The reason for choosing this
merit function is as follows: when the incident raychanges continuously in phase space along the phasespace boundary from point a to b, c, d, and back to a, thecorresponding ray at the image plane will also change
continuously in phase space and form a closed loop. According to the edge ray principle all the incident rays
within the rectangular region abcd in phase space will
fall into the closed loop in the phase space of the rays at
the image plane, bounded at the image plane within therange −ym to ym.
We choose the line aa⬘ to reach the ym point at the image plane, and also the line a⬘b to hit the upper edge ofba‘the lens back surface. Since a⬘ point has the smallest yep
value among those rays within the phase-space line a⬘b,
it should also have the largest ray angle when hitting theupper edge of the back surface so that its ray angle after
refraction by the back surface will also be largest. Thus, itwill reach the image plane with the largest y coordinate,ym, and all other rays along the line a⬘b in phase space
will reach the image plane below ym. Since the ray frompoint c in phase space hits the image plane at −ym, allrays from the phase-space line bc should reach the image
plane inside the concentration region 共−ym ⱕ y ⱕ ym兲. And,due to symmetry, the phase-space lines cd and da will fallinside the concentration region as well. Note that this approach is not exact: it is possible to generate a surface
which violates these assumptions. Moreover, while we canuse this one-dimensional design approach to generate a
two-dimensional (2D) surface by rotating the design surface about the axis, the analysis above has ignored skew
rays within the system. Nevertheless, we can obtain auseful design using the above approach, and upon optimizing M with the proximate ray trace equations for y
and y3, we obtain the design,

␣1 = 0.297 73,

␤1 = − 0.038 88,

␣2 = 0.018 92,

␤2 = 0.023 36,

␣3 = 0.001 71,

␤3 = 0.002 53,

Dep/2

共Dep/2兲−ya⬘+

Zheng et al.

+ ym

− ym

a
Fig. 6. (Color online) The approximate ray trace used by theedge ray principle for designing a concentrator lens. The labeledrays a, a⬘, and b are the phase-space points given in Fig. 5. Raya⬘ maps to the edge of the concentration region, at a distance ym
from the optical axis, while ray b is bent downward so that itreaches the image plane at y ⬍ ym.

The resulting lens is shown in Fig. 7 together with a diagram illustrating the ray mapping. Figure 8 shows the resulting concentration performance, giving the transmission (portion of rays reaching the concentration region) as
a function of the incidence angle.

8. CONCLUSION
The first implementation of proximate ray tracing, by
Hopkins in 1976 [6–8], appears to have been done as amethod of reducing the computational burden and complexity for calculating higher-order aberrations. By sampling a specific set of rays passing through the system,
one is able to obtain each of the various aberration coefficients. Our own implementation adapts the proximate ray
tracing concept for use in computer algebra systems in order to perform the entire procedure in the analytical domain. This is an important advance in that the analyticalformulas provide information of a qualitatively different
character than that provided by conventional ray tracing.
Some examples of these advantages include: (1) the aberration terms can be simply picked out of the final expression for the optical path length W = 兺snsws as a function of
the incidence angle 共Hx , Hy兲 and pupil coordinates
共xep , yep兲; (2) the optimization procedure can take advantage of the properties of well-behaved functions (e.g., polynomials) such as their infinite differentiability; (3) there
is no need to consider sampling density or similar issues
required for discrete ray tracing; and (4) the analytical
functional form more clearly shows some of the difficulties
that ray tracing can encounter, such as the presence of
singularities [25] that can affect the convergence of aber-Zheng et al.

Vol. 27, No. 8 / August 2010 / J. Opt. Soc. Am. A

1801

y
0.97
0.96
0.95
0.94
0.93

−1.5

−1.0

(a)

− 0.5

0.5

1.0

1.5

yep

(b)

Fig. 7. (Color online) (a) Layout of the concentrator lens obtained via the edge ray principle design method; (b) the profile of ray positiony at the image plane as a function of entrance pupil position yep for ␪ = 20°.

ration series. We plan to present a more thorough discussion of these properties in a subsequent publication.
Each of the advantages stated above can also be given
for the matrix approaches to analytical ray tracing[10–13]. The difference between the implementation presented here and the matrix methods is that our approach
is more easily adaptable to general-purpose use by having
the computer algebra program perform much of the generic analytical calculations, such as obtaining the surface normal, converting coordinate systems, and developing systems of equations out to arbitrary order. In
addition, while the research presented here considers
only monochromatic systems, we are currently adding
wavelength (or wavenumber, depending on the choice of
coordinate) as an additional primary variable and developing our code for use in multiwavelength systems. Thisis an essential additional step for optical system design
which will be new to analytical ray tracing.
The three design examples presented in Section 7 illustrate the flexibility of the analytical approach to treat systems of arbitrary symmetry and to produce accurate designs. In all three cases the optimization process in boththe proximate and numerical ray trace designs does not
Transmission
1.0

0.8

require interaction with the designer after defining thefirst-order properties of the system.
The drawbacks to analytical ray tracing include (1) forall but extreme cases, analytical ray tracing will be muchslower than numerical exact ray tracing; (2) the analytical formulas produced by the method can be quite
lengthy; and (3) analytical ray formulas are harder to interpret due to their unfamiliarity. While each of thesedrawbacks is important, they reflect the trade-off of the
new information obtained about a given design.

APPENDIX A
Due to the length of the polynomial expressions which
need to be manipulated for asymmetric systems and highorder approximations, an important feature of an analytical ray tracing engine is an efficient means of performing
each required operation. While [13] presents a code for
performing an order-expansion operation, the algorithm
used is quite slow for large expressions. An alternative
approach which is simple to apply for even very large expressions is the following. For each of the order-expanded
variables present in an expression, we can represent thesum as a vector where each element of the vector represents the appropriate order of expansion. Thus, for example, inside the computer algebra system we can represent the variable y as the vector

0.6

y = 共y共0兲,y共1兲,y共2兲,y共3兲, . . .兲T ,0.4

so that the multiplication operation between variables x
and y, for example, is computed via vector multiplications
as

0.2

共A1兲

Incident angles (degrees)
0

5

10

15

20

25

Fig. 8. (Color online) Concentrator example’s performance is illustrated by showing the portion of transmitted rays reaching
the concentration region as a function of incidence angle. The example shown here (solid line) compares well with that shown in
Fig. 8.9 of [24] (dashed line). A vertical line at 20° illustrates the
maximum angle of incidence used in the design.

1TxyT1,
where 1 is the vector of appropriate length containing 1’s
for each element. Here the inner operation, xyT, generates a matrix containing all combinations of elements of x
with elements of y, i.e.,

1802

J. Opt. Soc. Am. A / Vol. 27, No. 8 / August 2010

xyT =

冢

x 共0兲y 共0兲 x 共1兲y 共0兲 x 共2兲y 共0兲 . . .
x 共0兲y 共1兲 x 共1兲y 共1兲 x 共2兲y 共1兲 . . .
共0兲 共2兲

x y

共1兲 共2兲

x y

共2兲 共2兲

x y

...

]

]

]



冣

Zheng et al.
3.
4.

.

The outer operations perform the summation over all matrix elements. As given here, a computer algebra systemneeds to check the power of each of the variables used in
the expressions and add the powers of all variables withina given term, in order to obtain the order of a given term.
It is possible to achieve this in a single step by adding an
indicator variable, which we name g. The vector expression for g is given as

5.
6.
7.
8.
9.
10.

g = 共1,g,g ,g , . . .兲 .
2

3

T

When we define the vector representation of each variable
to be used, instead of Eq. (A1) we use

11.

y · g = 共y共0兲,gy共1兲,g2y共2兲,g3y共3兲, . . .兲,12.

where the “·” operator represents an element-by-element
multiplication. The indicator vector g is likewise applied
to all variables used in expressions applying orderexpansion methods. When operations such as the x-y multiplication shown above are performed, one now needs
only to locate the power of g in order to obtain the order ofa given term.
This ability is even more convenient for setting up andsolving systems of equations such as those for transfer or
refraction. For example, when calculating the expansionof the square root, as in Section 4, instead of manuallysetting up the matrix of equations and solving them at
each nonlinear operation, one needs only to ask the computer algebra system to perform the Taylor expansion of
the equation in the variable g. While not optimal in termsof the computational speed, this method greatly simplifies
the code, for improved readability and comprehension.

13.

REFERENCES
1.
2.

D. Shafer, “Global optimization in optical design,” Comput.
Phys. 8, 188–195 (1994).
A. E. W. Jones and G. W. Forbes, “An adaptive simulated
annealing algorithm for global optimization over continuous variables,” J. Global Optim. 6, 1–37 (1995).

14.
15.
16.

17.
18.
19.
20.
21.
22.
23.
24.
25.

S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, “Optimization by simulated annealing,” Science 220, 671–680 (1983).
X. Chen and K. Yamamoto, “An experiment in genetic optimization in lens design,” J. Mod. Opt. 44, 1693–1702(1997).
R. Horst, P. M. Pardalos, and N. V. Thoai, Introduction to
Global Optimization, 2nd ed. (Springer, 2000).
G. W. Hopkins, “Proximate ray tracing and optical aberration coefficients,” J. Opt. Soc. Am. 66, 405–410 (1976).
G. W. Hopkins, “Proximate ray tracing and wave aberration
coefficients,” J. Opt. Soc. Am. 66, 942–949 (1976).
G. W. Hopkins, “Aberrational analysis of optical systems: aproximate ray trace approach,” Ph.D. thesis (University of
Arizona, 1976).
J. von zur Gathen and J. Gerhard, Modern Computer Algebra (Cambridge U. Press, 1999).
M. Kondo and Y. Takeuchi, “Matrix method for nonlinear
transformation and its application to an optical lens system,” J. Opt. Soc. Am. A 13, 71–89 (1996).
V. Lakshminarayanan and S. Varadharajan, “Expressions
for aberrations coefficients using nonlinear transforms,”
Optom. Vision Sci. 74, 676–686 (1997).
J. B. Almeida, “General method for the determination of
matrix coefficients for high-order optical system modeling,”
J. Opt. Soc. Am. A 16, 596–601 (1999).
J. B. Almeida, “Programming matrix optics into Mathematica,” Optik (Stuttgart) 116, 270–276 (2005).
A. Walther, “Eikonal theory and computer algebra,” J. Opt.
Soc. Am. A 13, 523–531 (1996).
A. Walther, “Eikonal theory and computer algebra II,” J.
Opt. Soc. Am. A 13, 1763–1765 (1996).
T. Kryszczyński, “First steps towards an algebraic methodof the optical design in the range of all aberration orders,”
in 11th Slovak-Czech-Polish Optical Conference on Wave
and Quantum Aspects of Contemporary Optics (1999), pp.
336–342.
M. Born and E. Wolf, Principles of Optics, 7th ed. (Cambridge U. Press, 1999).
G. W. Forbes, “Optical system assessment for design: numerical ray tracing in the Gaussian pupil,” J. Opt. Soc. Am.
A 5, 1943–1956 (1988).
V. Y. Pan, “Solving a polynomial equation: some history andrecent progress,” SIAM Rev. 39, 187–220 (1997).
Wolfram Research, Inc., www.wolfram.com.
ZEMAX Development Corp., www.zemax.com.
D. J. Brady and N. Hagen, “Multiscale lens design,” Opt.
Express 17, 10,659–10,673 (2009).
N. Shatz and J. C. Bortz, “Global optimization of highperformance concentrators,” in Nonimaging Optics(Elsevier, 2005), Chap. 11, pp. 265–304.
R. Winston, J. C. Miñano, and P. Benítez, Nonimaging Optics (Elsevier, 2005).
G. W. Forbes, “Extension of the convergence of multivariateaberration series,” J. Opt. Soc. Am. A 3, 1376–1383 (1986).