Characterizing the shape of freeform optics

G. W. Forbes

doi:10.1364/OE.20.002483

1. Introduction

Progressive ophthalmic lenses were conceived over 75 years ago [1] and soon became more complex [2,3]. They clearly demonstrated the unique benefits that freeform optical surfaces can offer through the extra degrees of freedom available in surface shape. As documented in [4], other prominent commercial applications emerged in the 1970’s. On account of advances in optical fabrication and metrology, interest in a wider class of applications has escalated in recent decades. Within the domain of precision optics, these applications include head-mounted and head-up displays, phase elements for computational imaging and unobstructed reflective imaging systems (from ultra-compact digital cameras through projection television to EUV lithography). Lower-precision freeform surfaces have long been prominent in non-imaging applications such as illumination and it is impressive that some of the associated design methods have expanded in their capabilities so that they are now valuable for imaging applications as well [5,6]. Other interesting design methods are also evolving to enable the discovery of more effective systems involving freeform optics [7–9].

Standardized methods for communication of the nominal shape of freeforms are critical to this technology. Although no single convention can meet all the requirements of the various contexts and applications, some reasonably general-purpose options exist, e.g. point clouds with interpolation and various splines or wavelets [e.g. 10,11]. Alternatively, a common method in precision optics is to express the surface sag explicitly as a conic base with an additive deviation expressed as a polynomial. This has the advantage of being infinitely differentiable and is a simple generalization of the most widely used convention for characterizing rotationally symmetric optics. When the polynomial basis is orthogonalized, the associated coefficients become a spectrum of the surface’s frequency content that can be interpreted readily upon inspection. Further, such a spectrum can also facilitate estimates of manufacturability (which typically involve consideration of generating, testing, and oftentimes polishing). When implemented by using recurrence relations, this option is so robust and efficient that it enables an unlimited number of coefficients to be used as needed. These observations motivated the development of the polynomials introduced in this work.

Basic geometric ideas are presented in Section 2 where a general form is proposed to characterize the sag of a freeform optic. The polynomials used in this sag expression are derived in Section 3 and algorithms for their robust computation are presented in Appendices A and B. For demonstration, a particular surface is expressed in these terms in Section 4. A generalization is offered for certain special cases in Section 5 where some simple measures of testability are also discussed before the concluding remarks of Section 6.

2. Coordinates, aperture shapes, and sag expressions

Off-axis sections of rotationally symmetric surfaces can be regarded as a first step towards freeform optics. In such cases, it is natural to adopt coordinates that are aligned with the part’s axis of revolution. For example, the sag of a conventional conic section of revolution can then be expressed as

z = \frac{c_{con} [{(x + s)}^{2} + y^{2}]}{1 + \sqrt{1 - (1 + κ) c_{con}^{2} [{(x + s)}^{2} + y^{2}]}} = \frac{c_{con} (ρ^{2} + 2 s ρ \cos θ + s^{2})}{1 + \sqrt{1 - (1 + κ) c_{con}^{2} (ρ^{2} + 2 s ρ \cos θ + s^{2})}},

where

c_{con}

is the conic’s axial curvature and

κ

the conic constant. The first expression in Eq. (2.1) involves the Cartesian coordinates

{x, y, z}

with their origin displaced from the part’s vertex by

s

; the second is in terms of the associated cylindrical polar coordinates

{ρ, θ, z}

where

x = ρ \cos θ

and

y = ρ \sin θ

. Such a configuration with a vertical

z

axis is sketched in Fig. 1 .

Fig. 1 An off-axis section of a paraboloid is commonly encountered and has some of the challenges of more general freeform optics. The coordinate origin used for Eq. (2.1) is neither at the top nor the bottom of the cylinder as drawn for clarity at left, but in a plane that is tangent to the part’s vertex as shown at right. The separation between the cylinder and part axes is s.

Download Full Size | PDF

For freeform optics that are to be turned/ground, smoothed, polished, and tested —perhaps interferometrically with nulls or near-nulls and/or stitching— the more cost-effective shapes tend to be closer to spherical. Manufacturability typically grows in difficulty with departure from a best-fit sphere. More specifically, it is the rate of change of the departure along the normal from the best-fit sphere that must often be considered in both fabrication and testing. As discussed in [12], the difference between the local principal curvatures on the surface is often a driving consideration that can be related directly to derivatives of the normal departure. For many purposes, it is therefore helpful to characterize freeform shape by explicit reference to a best-fit sphere. That is one of the key motivators for the developments reported in what follows. Another of the critical considerations is that, for efficiency and numerical robustness, an orthogonal basis delivers significant benefits when characterizing shape [13,14].

It is not only the surface’s shape but also the shape of the part’s perimeter that must be characterized to fully specify a freeform. Rectangles, hexagons, ellipses, and kidney beans are familiar examples of nominal aperture shapes, but their specification is not always as simple as it may seem because the part’s perimeter typically does not lie in a plane. The aperture shown as the red curve in Fig. 1, for example, is defined as the intersection of a circular cylinder with the parent surface. For other geometries, the proposal in this work is to tightly enclose the part within a circular cylinder and then use the polar coordinates associated with that cylinder in the characterization of the surface’s shape, see Fig. 2 . Orthogonalization is necessarily performed over a specific domain and, for that purpose, I adopt this circular region that encloses the part.

Fig. 2 The characterization of a segment of a general freeform surface is depicted at left. In this case, as may be familiar from segmented mirrors, the segment of interest is taken to be the intersection of a hexagonal column with the surface. I opt for the coordinate origin to sit on the surface and choose the green cylinder to tightly enclose the segment of interest on the surface. The cylinder will typically be chosen to be nominally normal to the surface segment, as shown. The case for a part that resembles the paraboloidal segment of Fig. 1 is depicted at right.

Download Full Size | PDF

The next step therefore is to extend Eq. (2.1) of [14]. In particular, if we define $u$ by $u : = ρ / ρ_{\max}$ where $ρ_{\max}$ is the radius of the tightly enclosing cylinder, then the sag is now expressed as $z = f (ρ, θ)$ where

\begin{array}{l} f (ρ, θ) = \frac{c ρ^{2}}{1 + \sqrt{1 - c^{2} ρ^{2}}} + \frac{1}{\sqrt{1 - c^{2} ρ^{2}}} {u^{2} (1 - u^{2}) \sum_{n = 0}^{N} a_{n}^{0} Q_{n}^{0} (u^{2}) \\ + \sum_{m = 1}^{M} u^{m} \sum_{n = 0}^{N} [a_{n}^{m} \cos (m θ) + b_{n}^{m} \sin (m θ)] Q_{n}^{m} (u^{2})} . \end{array}

Here,

Q_{n}^{m} (v)

is a polynomial of order

n

in

v

,

c

is the curvature of what I somewhat loosely refer to as the “best-fit sphere”, and

\sqrt{1 - c^{2} ρ^{2}}

is the cosine of the angle between the

z

axis and the normal to that sphere. The presence of the cosine factor in Eq. (2.2) means that, to a first-order approximation, the entity within braces corresponds to a departure from the sphere that is measured along the local normals to that sphere.

Notice that the sum on the first line of Eq. (2.2) has precisely the same form as that used in [14] except that the coefficients and polynomials now carry a superscripted zero. This zero distinguishes them from the non-rotationally symmetric components, namely $a_{n}^{m}$ , $b_{n}^{m}$ , and $Q_{n}^{m} (u^{2})$ with $m > 0$ . Also note that the average over $θ$ of $f (ρ_{\max}, θ)$ is precisely the sag of the best-fit sphere at $ρ = ρ_{\max}$ : the component inside the braces averages to zero. In the unusual event of having to fit a predetermined shape with this representation, this observation can be useful to evaluate $c$ . More typically, all of the coefficients in Eq. (2.2) will be determined through system optimization during design.

Although the form of the polynomials is not locked down until Section 3, the similarity between the sum on the second line of Eq. (2.2) and the familiar decomposition into Zernike polynomials establishes that the terms in braces provide a complete set. As can be shown by using trigonometric identities, these terms are just a particular combination of the Cartesian elements $X^{j} Y^{k}$ where $X = u \cos θ$ and $Y = u \sin θ$ . [It is useful in this to equate real and imaginary parts after expanding the identity $u^{m} (\cos m θ + i \sin m θ) = {(X + i Y)}^{m}$ .] The maximum order of the elements that compose the contribution associated with $a_{n}^{m}$ or $b_{n}^{m}$ where $m > 0$ is just $m + 2 n$ , i.e. this is the maximum power of $u$ in those contributions. It would also be natural, therefore, to truncate the sums on the second line of Eq. (2.2) to include only the terms for which $m + 2 n \leq T$ for some $T$ . Due to the prefactor on the sum in the first line, where $m = 0$ , it is consistent for that sum to then extend up to $2 n + 4 \leq T$ .

Without loss of generality, it is possible to assume in Eq. (2.2) that

a_{0}^{1} = b_{0}^{1} = 0.

When they're non-zero, the contribution from these terms to the normal departure is proportional to

u (a_{0}^{1} \cos θ + b_{0}^{1} \sin θ) \propto (a_{0}^{1} x + b_{0}^{1} y)

. It is clear therefore that the cylinder's axis can be chosen to remove these average-tilt terms (much as the best-fit sphere was chosen to make the rotationally symmetric component vanish at the cylinder’s centre and edge). Alternatively, the cylinder can be oriented to remove either the mean tilt around its rim or the local tilt at the origin, which lead respectively to

\sum_{n = 0}^{N} a_{n}^{1} Q_{n}^{1} (1) = \sum_{n = 0}^{N} b_{n}^{1} Q_{n}^{1} (1) = 0, and \sum_{n = 0}^{N} a_{n}^{1} Q_{n}^{1} (0) = \sum_{n = 0}^{N} b_{n}^{1} Q_{n}^{1} (0) = 0.

Neither of Eqs. (2.4a) and (2.4b) is functionally as convenient as Eq. (2.3). I mention these options mainly to highlight some of the residual freedom in configuring the cylinder. Any of them can be used during optimization as part of design provided that the surface has appropriate tip and/or tilt parameters as separate degrees of freedom. In this way, it is assured that the base sphere of Eq. (2.2) is a decent approximation to the best-fit sphere as far as manufacturability is concerned. In the terminology of metrology, for example, the irrelevant piston, tilts, and power are suppressed in this way. The most important issue that remains is the process of constructing the new polynomials, i.e.

Q_{n}^{m} (v)

.

3. Orthogonalization in terms of the mean square gradient

Because it is the rate of change of the normal departure from a best-fit sphere that often drives manufacturability, the polynomials discussed in [12] were constructed to be orthogonal in slope. The mean square slope is then given by the sum of the squares of the coefficients in the departure. For freeform optics, it is the modulus of the gradient of the normal departure from the best-fit sphere that plays this role. That is, if the term in braces in Eq. (2.2) is written as

δ (u, θ) : = u^{2} (1 - u^{2}) \sum_{n = 0}^{N} a_{n}^{0} Q_{n}^{0} (u^{2}) + \sum_{m = 1}^{M} u^{m} \sum_{n = 0}^{N} [a_{n}^{m} \cos (m θ) + b_{n}^{m} \sin (m θ)] Q_{n}^{m} (u^{2}),

the natural extension of the basis used in [12] is to construct

Q_{n}^{m} (v)

so that the mean square gradient of this normal departure from the best-fit sphere satisfies

〈 {| \nabla δ (u, θ) |}^{2} 〉 = 〈 {(\frac{\partial δ}{\partial u})}^{2} + \frac{1}{u^{2}} {(\frac{\partial δ}{\partial θ})}^{2} 〉 = \sum_{m, n} [{(a_{n}^{m})}^{2} + {(b_{n}^{m})}^{2}] .

Here,

b_{n}^{0} = 0

and angle brackets denote the average over the aperture defined by

〈 g (u, θ) 〉 : = \int_{0}^{1} \int_{- π}^{π} g (u, θ) w (u) u d θ ​ d u / \int_{0}^{1} \int_{- π}^{π} w (u) u d θ ​ d u .

The weight function used in [12] is

w (u) = {[u^{2} (1 - u^{2})]}^{- 1 / 2}

so Eq. (3.3) becomes

〈 g (u, θ) 〉 : = \frac{1}{π^{2}} \int_{0}^{1} \int_{- π}^{π} g (u, θ) {[1 - u^{2}]}^{- 1 / 2} d θ ​ d u .

Equations (3.1), (3.2), and (3.4) completely determine

Q_{n}^{m} (v)

.

The trigonometric functions in Eq. (3.1) mean that the different azimuthal orders (i.e. $m θ$ and $m^{'} θ$ for $m \neq m^{'}$ ) are automatically orthogonal. That is, expanding the intermediate expression in Eq. (3.2) leads to sums of products, and the product of any two distinct azimuthal orders averages to zero. The only remaining step is to choose the polynomials within each separate azimuthal order so that they are internally orthonormal. Of course, $Q_{n}^{0} (v)$ are precisely the (“Qbfs”) polynomials introduced in [13], but those for $m > 0$ require further treatment. After performing the angle integral, the condition for orthonormality reduces to the requirement that

\frac{1}{π} \int_{0}^{1} {\frac{\partial}{\partial u} [u^{m} Q_{n}^{m} (u^{2})], m u^{m - 1} Q_{n}^{m} (u^{2})} \cdot {\frac{\partial}{\partial u} [u^{m} Q_{n^{'}}^{m} (u^{2})], m u^{m - 1} Q_{n^{'}}^{m} (u^{2})} \frac{d u}{\sqrt{1 - u^{2}}}

is unity when

n = n^{'}

and vanishes otherwise. With this, it is straightforward to construct the polynomials by using Gram-Schmidt orthogonalization. The results for

m

from 0 to 5 and

n

from 0 to 3 or 4 are presented in the tableau in Fig. 3 . Plots of the cosine versions of the associated functions in Eq. (3.1) along with contour plots of the modulus of their gradient are given in Fig. 4 . The complementary basis members, i.e.

u^{m} \sin (m θ) Q_{n}^{m} (u^{2})

, have similar plots, of course, but each is rotated by

π / (2 m)

.

Fig. 3 A sample of the polynomials constructed to be orthogonal in gradient. The columns for each azimuthal order contain the associated polynomials for n = 0, 1, 2,...

Download Full Size | PDF

Fig. 4 Plots of a sample of the basis members tabulated in Fig. 3. The four paired rows correspond to m = 0, 1, 2, and 5 while the five columns are for n = 0, 1, 2, 3, and 4. The upper row of each pair plots the associated basis function while the lower shows the modulus of their gradient. Importantly, the color scale is the same for all of the contour maps and it spans from 0 (purple) to 2 (red); a gradient of unity corresponds to the green of the plot for m = 1 and n = 0.

Download Full Size | PDF

The plots in Fig. 4 for $m = 0$ follow upon revolving the familiar plots in Fig. 1 of [14], although now it is the absolute value of the radial slope that is plotted. Recall that it is the non-uniform weight used in Eq. (3.4) that helps to constrain the peaks in the magnitude of the gradient to be not much bigger than $\sqrt{2}$ . This peak value is in keeping with a sine-like shape that has an rms value of unity. Also keep in mind that, as shown in Fig. 5 , the gradient is a vector and it is the integral of the inner product of these vector fields that underpins the orthonormalization. The polynomials listed in Fig. 3 allow for over 50 coefficients to be used in characterizing a freeform's shape [although the truncation discussed in the paragraph before Eq. (2.3) means that a subset of only $1 + 2 \times (3 + 2 + 2 + 1 + 1) = 19$ of these would be used if we stop at $T = 5$ , and Eqs. (2.3) and (2.4) may drop an additional two of these degrees of freedom regardless].

Fig. 5 Plots of the gradient field of three of the family of functions plotted in Fig. 4. It is the point-by-point inner products of different members of such maps that integrate to zero. The scale is fixed in all of these plots so that an arrow whose length matches the grid spacing represents a gradient of magnitude equal to the square root of two.

Download Full Size | PDF

To go beyond low orders, recurrence relations offer efficient, robust algorithms for evaluating arbitrary linear combinations of $Q_{n}^{m}$ and of their derivatives. One option to achieve this involves a link to Jacobi polynomials, and the associated details are presented in the Appendices. It is worth noting that a similar set of polynomials is developed in [15] to provide a basis for curl-free vector functions by way of a set of gradient-orthogonal polynomials. Interestingly, it is the different weight function [i.e. $w (u)$ in Eq. (3.3)] adopted here that causes the main difference with the polynomials derived in that work.

It is sometimes helpful to notice that the component in square brackets in Eq. (3.1) can trivially be re-expressed in terms of an alternative set of coefficients that give amplitude and phase:

a_{n}^{m} \cos (m θ) + b_{n}^{m} \sin (m θ) = α_{n}^{m} \cos (m θ - ϕ_{n}^{m}) .

With this, the mean square gradient of Eq. (3.2) satisfies

〈 {| \nabla δ (u, θ) |}^{2} 〉 = \sum_{m, n} [{(a_{n}^{m})}^{2} + {(b_{n}^{m})}^{2}] = \sum_{m, n} {(α_{n}^{m})}^{2} .

In this sense,

α_{n}^{m}

is then of primary importance in the spectrum while the phases become secondary. That is, a rough assessment of the dominant components of a part’s shape then follows upon inspection of just one half of the original number of coefficients.

4. Simple surface for demonstration

For demonstration and code verification, a simple paraboloid is expressed in the terms proposed above. This was also done in Section 4 of [14], but now an off-axis section of that same paraboloid is fitted with respect to a locally best-fit sphere. The paraboloid has curvature $1 / c_{parab} = 20 mm$ and the distance from its axis to its point of intersection with the centre line of the green cylinder drawn in Fig. 7 is taken to be $s = 20 mm$ . When the radius of the cylinder is $ρ_{\max} = 10 mm$ and the cylinder axis is chosen to be normal to the surface at its point of intersection, the curvature of the best-fit sphere is given by $1 / c = 37 .405283mm$ . (This curvature is given to an artificially high precision for code verification of a nm-level fit of the paraboloid with respect to this sphere.) The fit results are shown in Figs. 6 and 7 . Notice that symmetry in y means that $b_{n}^{m} \equiv 0$ . These results use some additional terms beyond those given in Fig. 3:

Fig. 6 Fitted coefficients in units of nm. These 23 coefficients are truncated as discussed in the paragraph before Eq. (2.3) with T = 8. If the part did not have a plane of symmetry, $b_{n}^{m}$ with m > 1 would bring the number of coefficients up to 43. These values have been rounded to the nearest nm and the results were used in computing the fit error of Fig. 7.

Download Full Size | PDF

Fig. 7 The sample paraboloid is shown at bottom along with the cylinder that encloses the off-axis section of interest and the associated best-fit sphere in red. The sag departure is plotted in the upper row along with the difference between this sag and that associated with the coefficients of Fig. 6. Notice that the fit error remains less than a nanometer.

Download Full Size | PDF

Q_{0}^{6} (x) = \frac{8}{9 \sqrt{7}}, Q_{1}^{6} (x) = \frac{8}{9 \sqrt{1397}} (77 - 72 x),

Q_{0}^{7} (x) = \frac{32}{7 \sqrt{231}}, Q_{0}^{8} (x) = \frac{8}{\sqrt{858}} .

While the coefficients in Fig. 6 correspond to a configuration satisfying Eq. (2.4b), the results for the configuration satisfying Eq. (2.3) are given in Fig. 8 . Notice the similarity between the two sets and that the shape is largely dominated by about $600 μ m$ of astigmatism ( $m = 2$ and $n = 0$ ) and $200 μ m$ of coma ( $m = 1$ and $n = 1$ ), see Fig. 4. That is, the character of the departure plot in Fig. 7 can be anticipated by a glance at these tables of coefficients.

Fig. 8 The corresponding coefficients for the case of Fig. 6, but where the cylinder has been tilted by 20.2223 mrad to enforce Eq. (2.3). That is, as highlighted above in red, $a_{0}^{1}$ now vanishes. The inverse of the best-fit curvature for this case is 37.432729 mm.

Download Full Size | PDF

5. The option for a non-zero conic constant

It is sometimes useful to express a part’s shape in terms of its deviation from a similar conicoid. This can be achieved by adding terms to Eq. (2.1) in a manner that generalizes Eq. (2.4) of [12]. In particular, the sag can be expressed as

z = \frac{c_{con} (ρ^{2} + 2 s ρ \cos θ + s^{2})}{1 + \sqrt{1 - (1 + κ) c_{con}^{2} (ρ^{2} + 2 s ρ \cos θ + s^{2})}} + \frac{δ (ρ / ρ_{\max}, θ)}{σ (ρ, θ)},

where the first term is just the conicoid’s sag of Eq. (2.1),

δ (u, θ)

was defined in Eq. (3.1), and

σ (ρ, θ)

is the cosine of the angle between the

z

axis and the normal to the conic, namely

σ (ρ, θ) : = \frac{\sqrt{1 - (1 + κ) c_{con}^{2} (ρ^{2} + 2 s ρ \cos θ + s^{2})}}{\sqrt{1 - κ c_{con}^{2} (ρ^{2} + 2 s ρ \cos θ + s^{2})}} .

To within a first-order approximation, the additive departure in sag from the conic in Eq. (5.1) can be converted to the departure along the normal to that conic by multiplying it by the cosine factor

σ (ρ, θ)

. It follows that

δ (u, θ)

itself can again be regarded as the departure along the normal to the reference surface. This generalization enables the handling of extremely fast cases that would otherwise involve unworkable hyper-hemispheres. All the methods derived in the appendices carry over without change because, once again, it is the sum in Eq. (B.2) that stands as the core component. This is also an effective option when characterizing parts that are close in shape to a conic section for which an interferometric null test is in hand and where the metrology represents a critical part of manufacturability.

A part that approximates a paraboloid, for example, can be tested by using a retro sphere and a transmission flat as depicted in Fig. 1 of [12], but now imagine that only the upper half of that figure holds the off-axis segment of interest. As described in [12], a map of the fringe density in the resulting interferogram is approximately proportional to $| \nabla δ (u, θ) |$ times the cosine of the angle of incidence of the test rays at the part. In many cases, the cosine factor is roughly constant. While it is the peak fringe density that ultimately drives testability, the rms fringe density is a useful and more readily accessible measure. Thankfully, the weight function of Eq. (3.4) serves to constrain the ratio of peak to rms values. When such a test is performed at wavelength $λ$ on an $N \times N$ pixel grid, it follows from the discussion in [12] that the rms fringe density in Nyquist units, say $\bar{γ}$ , satisfies

\bar{γ} \approx \frac{16}{N λ} \sqrt{\sum_{m, n} [{(a_{n}^{m})}^{2} + {(b_{n}^{m})}^{2}]} .

Results such as this can be used to derive simple design constraints and valuable estimates of testability.

It is striking that Eq. (5.3) involves little more than summing the squares of the shape coefficients; there is no need to evaluate a single departure from a best-fit conic or to average its squared rate of change. The orthogonalization introduced in Section 3 was constructed to make this rms gradient trivially accessible. A similar result applies, of course, for a conventional non-null interferometric test with a transmission sphere for a part described by using Eq. (2.2). [The single bounce off the part in that case means that the factor of 16 in Eq. (5.3) is replaced by an 8.] What is more, the same analysis also applies when considering the line spacing of a CGH null. Do not forget, however, that the sum of the squares in Eq. (5.3) gives an rms over the enclosing cylinder rather than just the part’s clear aperture. In the later stages of design or fabrication, more accurate analyses will typically be required.

6. Concluding remarks

The diversity of applications of freeform optics means that no single convention for characterizing shape will be ideal for all cases. For example, when a part has a rectangular aperture and is manufactured by rastering processes, it would be natural to construct a description in terms of Cartesian coordinates rather than the cylindrical polar coordinates used above. Alternatively, for parts that are ready to be deployed immediately after generation by single-point diamond turning, it is the azimuthal component of the gradient of the sag that oftentimes drives manufacturability, and direct access to this property could be factored into their characterization. However, for high-precision optics that must be turned/ground, smoothed, finished, and perhaps interferometrically tested (probably with nulls or near-nulls and/or stitching), it is effective to express their shape with respect to a best-fit sphere. In particular, it is typically the gradient of the normal departure from the best-fit sphere that is of crucial importance. This observation motivated the generalization of Qbfs (of [13]) that was presented in Sections 2 and 3 in terms of an enclosing cylinder. The further generalizations offered in Section 5 include not only a conic constant, but also the option to move the base surface’s center of curvature away from the cylinder’s axis. This second freedom can be exploited even when . For general purposes, however, it is Eq. (2.2) that is preferred over Eq. (5.1) because other aspects of manufacturability are more readily assessed when a freeform surface is expressed with respect to such a best-fit sphere.

When supplemented with the powerful algorithms presented in the appendices, Eq. (2.2) offers a general-purpose option for characterizing freeform shape in terms of coefficients that form an intuitive and elegant spectrum. For ray tracing through such a freeform, it can be helpful to have explicit access to the surface’s local normal vector. In terms of coordinates related by ${x, y, z} \equiv {ρ \cos θ, ρ \sin θ, z}$ , the Cartesian components of the normal vector for the surface described by $z = f (ρ, θ)$ are readily found to be

{f_{ρ} \cos θ - \frac{1}{ρ} f_{θ} \sin θ, f_{ρ} \sin θ + \frac{1}{ρ} f_{θ} \cos θ, - 1},

where subscripts denote partial derivatives. Similarly, for first-order aberration analysis, it can be helpful to note that this surface can be approximated near the coordinate origin by

\begin{array}{l} z = [\frac{1}{ρ_{\max}} \sum_{n} a_{n}^{1} Q_{n}^{1} (0)] x + [\frac{1}{ρ_{\max}} \sum_{n} b_{n}^{1} Q_{n}^{1} (0)] y + [\frac{1}{2} c + \frac{1}{ρ_{\max}^{2}} \sum_{n} a_{n}^{0} Q_{n}^{0} (0)] (x^{2} + y^{2}) \\ + [\frac{1}{ρ_{\max}^{2}} \sum_{n} a_{n}^{2} Q_{n}^{2} (0)] (x^{2} - y^{2}) + 2 [\frac{1}{ρ_{\max}^{2}} \sum_{n} b_{n}^{2} Q_{n}^{2} (0)] x y + O (3) . \end{array}

The results presented in the appendices give simple access to these linear combinations of

Q_{n}^{m} (0)

for fixed

m > 0

, see Eqs. (A.23) and (B.2) to (B.4). Notice that this approach is a direct analog of Eq. (3.16) of [14], which can now be re-used for the sum with

m = 0

in Eq. (6.2). This further highlights the important observation that Eq. (2.2) includes as a special case the option for characterizing rotationally symmetric aspheres that was introduced in [13]. Because it is built on notions of manufacturability, this single framework can therefore meet many of both the traditional as well as the emerging needs of our industry.

Appendix A: Auxiliary polynomials

Be warned at the outset that the equations in this opening paragraph will be modified below in order to avoid an isolated and surprising failure. The key here is that the polynomials defined in Section 3 can be expressed in terms of the Jacobi polynomials. In particular, consider the auxiliary polynomials defined by

P_{n}^{m} (x) : = {(- 1)}^{n} \frac{(2 n)!!}{2 (2 n - 1)!!} P_{n}^{(- \frac{3}{2}, m - \frac{3}{2})} (2 x - 1),

where $P_{n}^{(α, β)} (v)$ is the Jacobi polynomial of order $n$ . Here, the factorial operators are defined by $0! = 1$ with $n! = n (n - 1)!$ and $(- 1)!! = 0!! = 1$ with $n!! = n (n - 2)!!$ . The polynomials of Eq. (A.1) satisfy a standard three-term recurrence relation of the form

P_{n + 1}^{m} (x) = [A_{n}^{m} + B_{n}^{m} x] P_{n}^{m} (x) - C_{n}^{m} P_{n - 1}^{m} (x),

where the constants are given explicitly by

A_{n}^{m} : = (2 n - 1) (m + 2 n - 2) [4 n (m + n - 2) + (m - 3) (2 m - 1)] / D_{n}^{m},

B_{n}^{m} : = - 2 (2 n - 1) (m + 2 n - 3) (m + 2 n - 2) (m + 2 n - 1) / D_{n}^{m},

C_{n}^{m} : = n (2 n - 3) (m + 2 n - 1) (2 m + 2 n - 3) / D_{n}^{m},

with

D_{n}^{m} : = (4 n^{2} - 1) (m + n - 2) (m + 2 n - 3) .

Just as the auxiliary basis did for the polynomials used in [14], $P_{n}^{m} (x)$ and Eq. (A.2) underpin a number of efficient processes for working with $Q_{n}^{m} (x)$ . The scaling factor in Eq. (A.1) is chosen so that the peak absolute value of $u^{m} P_{n}^{m} (u^{2})$ over $0 \leq u \leq 1$ never exceeds about unity for all m and n. As in [14], this makes it straightforward to determine when terms are negligible and can be dropped from a linear combination.

For any $m \geq 1$ , the recurrence relation in Eq. (A.2) would normally be initialized by using $P_{0}^{m} (x) = 1 / 2$ and $P_{1}^{m} (x) = m - \frac{1}{2} - (m - 1) x$ . [In fact, just the first of these is sufficient for $m > 3$ since $C_{0}^{m} = 0$ and Eq. (A.2) itself generates $P_{1}^{m} (x)$ by taking $n = 0$ .] Notice that for $m = 1$ , however, $P_{1}^{m} (x)$ would then not be a linear function but a constant. Further, $D_{n}^{m}$ vanishes when $m = n = 1$ thereby causing Eqs. (A.2) and (A.3) to fail in the generation of $P_{2}^{1} (x)$ . Special handling is therefore required and one option is just to patch in a special clause. In what follows, I therefore replace Eq. (A.1) by

P_{n}^{m} (x) : = {\begin{cases} 1 - x / 2, m = 1 and n = 1, \\ \frac{{(- 1)}^{n} (2 n)!!}{2 (2 n - 1)!!} P_{n}^{(- \frac{3}{2}, m - \frac{3}{2})} (2 x - 1), otherwise . \end{cases}

In this case, the first two polynomials are

P_{0}^{m} (x) = \frac{1}{2} and P_{1}^{m} (x) = {\begin{cases} 1 - x / 2, m = 1, \\ m - \frac{1}{2} - (m - 1) x, m > 1. \end{cases}

This patch means that, when $m = 1$ , Eq. (A.2) is valid only for $n > 2$ , and the initialization for that case is then

P_{2}^{1} (x) = [3 - x (12 - 8 x)] / 6 and P_{3}^{1} (x) = {5 - x [60 - x (120 - 64 x)]} / 10.

If the integral in Eq. (3.5) is taken to define the inner product written as $〈 Q_{n}^{m}, Q_{n^{'}}^{m} 〉$ then it turns out that, for any fixed $m \geq 1$ , $〈 P_{n}^{m}, P_{n^{'}}^{m} 〉$ is a symmetric tri-diagonal matrix. The proof of this follows upon applying five standard identities, namely

\frac{\partial}{\partial x} [P_{n}^{(a, b)} (2 x - 1)] = (n + a + b + 1) P_{n - 1}^{(a + 1, b + 1)} (2 x - 1),

P_{n}^{(- \frac{3}{2}, m - \frac{3}{2})} (z) = \frac{1}{m + 2 n - 2} [(m + n - 2) P_{n}^{(- \frac{1}{2}, m - \frac{3}{2})} (z) - (m + n - \frac{3}{2}) P_{n - 1}^{(- \frac{1}{2}, m - \frac{3}{2})} (z)],

\int_{0}^{1} x^{m - \frac{3}{2}} {[P_{n}^{(- \frac{1}{2}, m - \frac{3}{2})} (2 x - 1)]}^{2} \frac{d x}{\sqrt{1 - x}} = π \frac{(2 n - 1)!! (2 m + 2 n - 3)!!}{2^{m + 2 n - 1} (m + 2 n - 1) n! (m + n - 2)!},

P_{n}^{(- \frac{1}{2}, m - \frac{3}{2})} (z) = \frac{1}{m + 2 n - 1} [(m + n - 1) P_{n}^{(- \frac{1}{2}, m - \frac{1}{2})} (z) + (n - \frac{1}{2}) P_{n - 1}^{(- \frac{1}{2}, m - \frac{1}{2})} (z)],

\begin{array}{l} x P_{n}^{(- \frac{1}{2}, m - \frac{1}{2})} (2 x - 1) = {\frac{(n + 1) (m + n)}{(m + 2 n + 1) (m + 2 n)} P_{n + 1}^{(- \frac{1}{2}, m - \frac{1}{2})} (2 x - 1) \\ + \frac{[4 m n + m (2 m - 1) + 4 n^{2} - 1]}{2 (m + 2 n + 1) (m + 2 n - 1)} P_{n}^{(- \frac{1}{2}, m - \frac{1}{2})} (2 x - 1) \\ + \frac{(2 n - 1) (2 m + 2 n - 1)}{4 (m + 2 n) (m + 2 n - 1)} P_{n - 1}^{(- \frac{1}{2}, m - \frac{1}{2})} (2 x - 1)} . \end{array}

This process involves replacing whichever of the expressions on the left-hand sides of Eqs. (A.8)–(A.12) appears first by the associated right-hand side, and the final step is to invoke Eq. (A.10) once more but with $m$ replaced throughout by $m + 1$ . While this sounds somewhat daunting, present-day computer algebra is up to the task when provided some guidance. In this way, it follows from the orthogonality of the Jacobi polynomials that $〈 P_{n}^{m}, P_{n^{'}}^{m} 〉$ is zero whenever $| n - n^{'} | > 1$ . The diagonal elements, say $F_{n}^{m} : = 〈 P_{n}^{m}, P_{n}^{m} 〉$ , turn out to be given by

F_{n}^{m} = {\begin{cases} \frac{m^{2} (2 m - 3)!!}{2^{m + 1} (m - 1)!}, n = 0, \\ \frac{4 {(n - 1)}^{2} n^{2} + 1}{8 {(2 n - 1)}^{2}} + \frac{11}{32} δ_{n 1}, n > 0 and m = 1, \\ \frac{2 n χ (3 - 5 m + 4 n χ) + m^{2} (3 - m + 4 n χ)}{(m + 2 n - 3) (m + 2 n - 2) (m + 2 n - 1) (2 n - 1)} γ_{n}^{m}, n > 0 and m > 1, \end{cases}

where $χ = m + n - 2$ as a shorthand, $δ_{m n}$ is the Kronecker delta, and

γ_{n}^{m} : = \frac{n! (2 m + 2 n - 3)!!}{2^{m + 1} (m + n - 3)! (2 n - 1)!!} .

Similarly, the off-diagonal elements, say $G_{n}^{m} : = 〈 P_{n}^{m}, P_{n + 1}^{m} 〉$ , are found to be

G_{n}^{m} = {\begin{cases} \frac{(2 m - 1)!!}{2^{m + 1} (m - 1)!}, n = 0, \\ - \frac{(2 n^{2} - 1) (n^{2} - 1)}{8 (4 n^{2} - 1)} - \frac{1}{24} δ_{n 1}, n > 0 and m = 1, \\ - \frac{[2 n (m + n - 1) - m] (n + 1) (2 m + 2 n - 1)}{(m + 2 n - 2) (m + 2 n - 1) (m + 2 n) (2 n + 1)} γ_{n}^{m}, n > 0 and m > 1. \end{cases}

Notice that $γ_{n}^{m}$ of Eq. (A.14) is required only for $n > 0$ and $m > 1$ , and this can be computed effectively by starting with $γ_{1}^{2} = 3 / 8$ and working up progressively, first in m and then in n, by using

γ_{1}^{m + 1} = \frac{2 m + 1}{2 (m - 1)} γ_{1}^{m} and γ_{n + 1}^{m} = \frac{(n + 1) (2 m + 2 n - 1)}{(m + n - 2) (2 n + 1)} γ_{n}^{m} .

Of course, a similar process can be used for evaluating $F_{0}^{m}$ and $G_{0}^{m}$ .

The simple and powerful consequence of the tri-diagonal result derived in the previous paragraph is that —much as described at Eq. (A.12) of [14]— $Q_{n}^{m} (x)$ becomes accessible through a recurrence relation:

P_{n}^{m} (x) = {\begin{cases} f_{0}^{m} Q_{0}^{m} (x), n = 0, \\ f_{n}^{m} Q_{n}^{m} (x) + g_{n - 1}^{m} Q_{n - 1}^{m} (x), n > 0, \end{cases}

where the constants $f_{n}^{m}$ and $g_{n}^{m}$ are found by Cholesky decomposition of $〈 P_{n}^{m}, P_{n^{'}}^{m} 〉$ . That is, for each azimuthal order —i.e. for any fixed m— these constants are determined by starting with $f_{0}^{m} = \sqrt{F_{0}^{m}}$ and working progressively up from $n = 1$ to any desired terminal value by using

g_{n - 1}^{m} = G_{n - 1}^{m} / f_{n - 1}^{m}, and f_{n}^{m} = \sqrt{F_{n}^{m} - g_{n - 1}^{m} g_{n - 1}^{m}} .

The significance of Eqs. (A.2) and (A.17) is discussed further in Appendix B. For efficiency, $A_{n}^{m}$ , $B_{n}^{m}$ , $C_{n}^{m}$ , $f_{n}^{m}$ and $g_{n}^{m}$ can easily be pre-computed and stored over any desired range of m and n. For verification of any implementation of these processes, I note that the initial sub-blocks of all the entities appearing in Eq. (A.18) are given by

[\begin{matrix} F_{0}^{1} \\ F_{1}^{1} \\ F_{2}^{1} \\ F_{3}^{1} \end{matrix} \begin{matrix} F_{0}^{2} \\ F_{1}^{2} \\ F_{2}^{2} \\ F_{3}^{2} \end{matrix} \begin{matrix} F_{0}^{3} \\ F_{1}^{3} \\ F_{2}^{3} \\ F_{3}^{3} \end{matrix} \begin{matrix} F_{0}^{4} \\ F_{1}^{4} \\ F_{2}^{4} \\ F_{3}^{4} \end{matrix}] = [\begin{matrix} \frac{1}{4} \\ \frac{15}{32} \\ \frac{17}{72} \\ \frac{29}{40} \end{matrix} \begin{matrix} \frac{1}{2} \\ \frac{7}{8} \\ \frac{35}{36} \\ \frac{67}{40} \end{matrix} \begin{matrix} \frac{27}{32} \\ \frac{35}{16} \\ \frac{35}{16} \\ \frac{243}{80} \end{matrix} \begin{matrix} \frac{5}{4} \\ \frac{511}{128} \\ \frac{23}{6} \\ \frac{12287}{2560} \end{matrix}], [\begin{matrix} G_{0}^{1} \\ G_{1}^{1} \\ G_{2}^{1} \end{matrix} \begin{matrix} G_{0}^{1} \\ G_{1}^{2} \\ G_{2}^{2} \end{matrix} \begin{matrix} G_{0}^{1} \\ G_{1}^{3} \\ G_{2}^{3} \end{matrix} \begin{matrix} G_{0}^{1} \\ G_{1}^{4} \\ G_{2}^{4} \end{matrix}] = [\begin{matrix} \frac{1}{4} \\ \frac{1}{24} \\ \frac{7}{40} \end{matrix} \begin{matrix} \frac{3}{8} \\ \frac{5}{48} \\ \frac{7}{16} \end{matrix} \begin{matrix} \frac{15}{32} \\ \frac{7}{32} \\ \frac{117}{160} \end{matrix} \begin{matrix} \frac{35}{64} \\ \frac{21}{64} \\ \frac{33}{32} \end{matrix}],

[\begin{matrix} f_{0}^{1} & f_{0}^{2} & f_{0}^{3} & f_{0}^{4} \\ f_{1}^{1} & f_{1}^{2} & f_{1}^{3} & f_{1}^{4} \\ f_{2}^{1} & f_{2}^{2} & f_{2}^{3} & f_{2}^{4} \\ f_{3}^{1} & f_{3}^{2} & f_{3}^{3} & f_{3}^{4} \end{matrix}] = [\begin{matrix} \frac{1}{2} & \frac{1}{\sqrt{2}} & \frac{3}{4} \sqrt{\frac{3}{2}} & \frac{\sqrt{5}}{2} \\ \frac{1}{4} \sqrt{\frac{7}{2}} & \frac{1}{4} \sqrt{\frac{19}{2}} & \frac{1}{4} \sqrt{\frac{185}{6}} & \frac{3 \sqrt{427}}{32} \\ \frac{1}{6} \sqrt{\frac{115}{14}} & \frac{1}{2} \sqrt{\frac{145}{38}} & \frac{1}{4} \sqrt{\frac{12803}{370}} & \frac{1}{2} \sqrt{\frac{2785}{183}} \\ \frac{1}{5} \sqrt{\frac{3397}{230}} & \frac{1}{4} \sqrt{\frac{6841}{290}} & \frac{9}{4} \sqrt{\frac{14113}{25606}} & \frac{1}{16} \sqrt{\frac{1289057}{1114}} \end{matrix}],

[\begin{matrix} g_{0}^{1} & g_{0}^{2} & g_{0}^{3} & g_{0}^{4} \\ g_{1}^{1} & g_{1}^{2} & g_{1}^{3} & g_{1}^{4} \\ g_{2}^{1} & g_{2}^{2} & g_{2}^{3} & g_{2}^{4} \end{matrix}] = [\begin{matrix} \frac{1}{2} & \frac{3}{4 \sqrt{2}} & \frac{5}{4 \sqrt{6}} & \frac{7 \sqrt{5}}{32} \\ - \frac{1}{3 \sqrt{14}} & - \frac{5}{6 \sqrt{38}} & - \frac{7}{4} \sqrt{\frac{3}{370}} & - \frac{1}{2} \sqrt{\frac{7}{61}} \\ - \frac{21}{10} \sqrt{\frac{7}{230}} & - \frac{7}{4} \sqrt{\frac{19}{290}} & - \frac{117}{4} \sqrt{\frac{37}{128030}} & - \frac{33}{16} \sqrt{\frac{183}{2785}} \end{matrix}] .

The results in this appendix give ready access to $Q_{n}^{m} (x)$ at any order. In particular, it follows from Eq. (A.17) that we can start with $Q_{0}^{m} (x) = P_{0}^{m} (x) / f_{0}^{m} = 1 / (2 f_{0}^{m})$ and then work up in $n$ by using

Q_{n}^{m} (x) = [P_{n}^{m} (x) - g_{n - 1}^{m} Q_{n - 1}^{m} (x)] / f_{n}^{m} .

Keep in mind that $P_{n}^{m} (x)$ is accessible through Eqs. (A.2), (A.3), (A.5), and (A.6). More interestingly, however, it is shown in Appendix B that a linear combination of either of these polynomial sets can be evaluated directly without computing any of the polynomials themselves. In fact, the methods in Appendix B also offer cleaner ways to evaluate $Q_{n}^{m} (x)$ . It is also helpful to note for the purposes of first-order aberration analysis in cases where the chief ray intersects the freeform at the axis of the cylinder, that it follows from Eq. (A.4) that

P_{n}^{m} (0) = {\begin{cases} 1, m = 1 and n = 1, \\ \frac{(2 m + 2 n - 3)!!}{2 (2 m - 3)!! (2 n - 1)!!}, otherwise . \end{cases}

Much as Eq. (A.1) was ultimately replaced by Eq. (A.4), the warning at the outset of this appendix carries into the next one, which involves modifications of Eqs. (A.3). Keep in mind that these admittedly ungainly patches have no impact on $Q_{n}^{m} (x)$ at all; they enter only the (typically hidden) auxiliary computational processes offered in these two appendices.

Appendix B: Additional recurrence-based processes

When working with sag functions defined by Eq. (2.2) it is helpful to re-express this result as

\begin{array}{l} f (ρ, θ) = \frac{c ρ^{2}}{1 + \sqrt{1 - c^{2} ρ^{2}}} + \frac{1}{\sqrt{1 - c^{2} ρ^{2}}} {u^{2} (1 - u^{2}) \sum_{n = 0}^{N} a_{n}^{0} Q_{n}^{0} (u^{2}) \\ + \sum_{m = 1}^{M} u^{m} [\cos (m θ) \sum_{n = 0}^{N} a_{n}^{m} Q_{n}^{m} (u^{2}) + \sin (m θ) \sum_{n = 0}^{N} b_{n}^{m} Q_{n}^{m} (u^{2})]} . \end{array}

Efficient processes to evaluate the sum on the first line of Eq. (B.1) and to determine its derivatives of any order were presented in [14]. For any fixed $m \geq 1$ , this appendix addresses the additional core task of working with the sums of the form of those on the second line, say

S (x) : = \sum_{n = 0}^{N} c_{n} Q_{n}^{m} (x),

where $c_{n}$ may be either $a_{n}^{m}$ or $b_{n}^{m}$ . (To reduce clutter, neither $S$ nor $c_{n}$ takes a superscripted $m$ .) By following the same procedure discussed in Section 3 of [14], it follows from Eq. (A.17) above that $S (x)$ can also be expressed in terms of the auxiliary polynomials as

S (x) \equiv \sum_{n = 0}^{N} d_{n} P_{n}^{m} (x),

where $d_{n}$ can be found simply from $c_{n}$ by starting with $d_{N} = c_{N} / f_{N}^{m}$ and working down in $n$ from $N - 1$ to 0 with

d_{n} = (c_{n} - g_{n}^{m} d_{n + 1}) / f_{n}^{m} .

The inverse of this process is also sometimes useful: for $n$ from 0 to $N - 1$ use

c_{n} = f_{n}^{m} d_{n} + g_{n}^{m} d_{n + 1}

and finish with $c_{N} = f_{N}^{m} d_{N}$ . That is, for any $m > 0$ , superpositions of $Q_{n}^{m} (x)$ and $P_{n}^{m} (x)$ are trivially interchangeable in this way.

To avoid complications associated with the patch in Eq. (A.4), first consider working with Eq. (B.3) for any $m > 1$ . For these cases, the recurrence relation in Eq. (A.2) means that the processes discussed in [16] are directly applicable. In particular, $S (x)$ can be evaluated by using Clenshaw’s method, which starts with $α_{N} = d_{N}$ and $α_{N - 1} = d_{N - 1} + (A_{N - 1}^{m} + B_{N - 1}^{m} x) α_{N}$ and progressively works down from $n = N - 2$ with

α_{n} = d_{n} + (A_{n}^{m} + B_{n}^{m} x) α_{n + 1} - C_{n + 1}^{m} α_{n + 2} .

(An equivalent initialization is $α_{N + 2} = α_{N + 1} = 0$ and now work down from $n = N$ .) The desired end result is given simply by $S (x) = \frac{1}{2} α_{0}$ .

When $m = 1$ , on the other hand, it helps to choose values for $A_{n}^{1}$ , $B_{n}^{1}$ , and $C_{n}^{1}$ , for $n$ equal to 0, 1, and 2 so that Eq. (A.2) is satisfied. This can only partially be achieved for $n = 2$ where just three of the four coefficients in the cubics on the left- and right-hand sides can be matched. The discrepancy is chosen to be a constant here, and it turns out to be $- 2 / 5$ . With this, when Eqs. (A.3) are modified to accept the following special cases

A_{0}^{1} = 2, B_{0}^{1} = - 1, A_{1}^{1} = - 4 / 3, B_{1}^{1} = - 8 / 3, C_{1}^{1} = - 11 / 3, C_{2}^{1} = 0,

A_{0}^{m} = 2 m - 1, B_{0}^{m} = 2 (1 - m), when m > 1,

it can be shown that the desired result is then given in all cases by

S (x) = {\begin{cases} \frac{1}{2} α_{0} - \frac{2}{5} α_{3}, m = 1 and N > 2, \\ \frac{1}{2} α_{0} otherwise . \end{cases}

Notice that $α_{n}$ of Eq. (B.6) is a polynomial in $x$ of order $N - n$ .

One of the most striking aspects of Eqs. (B.6) and (B.9) is that, as described in Appendix B of [16], they allow derivatives of any order to be evaluated with similar ease. If $S^{(j)} (x)$ denotes the $j$ ’th derivative of $S (x)$ then it follows from Eqs. (B.9) and (B.6) respectively that

S^{(j)} (x) = {\begin{cases} \frac{1}{2} α_{0}^{(j)} - \frac{2}{5} α_{3}^{(j)}, m = 1 and N > 2, \\ \frac{1}{2} α_{0}^{(j)} otherwise, \end{cases}

α_{n}^{(j)} = j B_{n}^{m} α_{n + 1}^{(j - 1)} + (A_{n}^{m} + B_{n}^{m} x) α_{n + 1}^{(j)} - C_{n + 1}^{m} α_{n + 2}^{(j)}, when j > 0.

When initialized with $α_{N + 2 - j}^{(j)} = α_{N + 1 - j}^{(j)} = 0$ , Eq. (B.11) can proceed downward from $n = N - j$ to give the desired results for Eq. (B.10). Notice that this process for evaluating the $j$ ’th derivative assumes that the $(j - 1)$ ’th derivative has already been determined. Also, it is often the case that a particular $S (x)$ of Eq. (B.2) may need to be evaluated at many different values of $x$ . It is worth noting therefore that the conversion from $c_{n}$ to $d_{n}$ of Eq. (B.4) need then be done only once. It is Eqs. (B.6), (B.10), and (B.11) that give efficient access to $S (x)$ and all of its derivatives. Finally, for code verification, I note that the low-order blocks of the associated entities after the patches of Eqs. (B.7) and (B.8) (shown below in red) are

[\begin{matrix} A_{0}^{1} & A_{0}^{2} & A_{0}^{3} & A_{0}^{4} & A_{0}^{5} & A_{0}^{6} \\ A_{1}^{1} & A_{1}^{2} & A_{1}^{3} & A_{1}^{4} & A_{1}^{5} & A_{1}^{6} \\ A_{2}^{1} & A_{2}^{2} & A_{2}^{3} & A_{2}^{4} & A_{2}^{5} & A_{2}^{6} \\ A_{3}^{1} & A_{3}^{2} & A_{3}^{3} & A_{3}^{4} & A_{3}^{5} & A_{3}^{6} \\ A_{4}^{1} & A_{4}^{2} & A_{4}^{3} & A_{4}^{4} & A_{4}^{5} & A_{4}^{6} \end{matrix}] = [\begin{matrix} 2 & 3 & 5 & 7 & 9 & 11 \\ - \frac{4}{3} & \frac{2}{3} & 2 & \frac{76}{27} & \frac{85}{24} & \frac{106}{25} \\ \frac{9}{5} & \frac{26}{15} & 2 & \frac{117}{50} & \frac{203}{75} & \frac{108}{35} \\ \frac{55}{28} & \frac{66}{35} & 2 & \frac{536}{245} & \frac{135}{56} & \frac{130}{49} \\ \frac{161}{81} & \frac{122}{63} & 2 & \frac{515}{243} & \frac{143}{63} & \frac{161}{66} \end{matrix}],

[\begin{matrix} B_{0}^{1} & B_{0}^{2} & B_{0}^{3} & B_{0}^{4} & B_{0}^{5} & B_{0}^{6} \\ B_{1}^{1} & B_{1}^{2} & B_{1}^{3} & B_{1}^{4} & B_{1}^{5} & B_{1}^{6} \\ B_{2}^{1} & B_{2}^{2} & B_{2}^{3} & B_{2}^{4} & B_{2}^{5} & B_{2}^{6} \\ B_{3}^{1} & B_{3}^{2} & B_{3}^{3} & B_{3}^{4} & B_{3}^{5} & B_{3}^{6} \\ B_{4}^{1} & B_{4}^{2} & B_{4}^{3} & B_{4}^{4} & B_{4}^{5} & B_{4}^{6} \end{matrix}] = [\begin{matrix} - 1 & - 2 & - 4 & - 6 & - 8 & - 10 \\ - \frac{8}{3} & - 4 & - 4 & - \frac{40}{9} & - 5 & - \frac{28}{5} \\ - \frac{24}{5} & - 4 & - 4 & - \frac{21}{5} & - \frac{112}{25} & - \frac{24}{5} \\ - \frac{30}{7} & - 4 & - 4 & - \frac{144}{35} & - \frac{30}{7} & - \frac{220}{49} \\ - \frac{112}{27} & - 4 & - 4 & - \frac{110}{27} & - \frac{88}{21} & - \frac{13}{3} \end{matrix}],

[\begin{matrix} C_{0}^{1} & C_{0}^{2} & C_{0}^{3} & C_{0}^{4} & C_{0}^{5} & C_{0}^{6} \\ C_{1}^{1} & C_{1}^{2} & C_{1}^{3} & C_{1}^{4} & C_{1}^{5} & C_{1}^{6} \\ C_{2}^{1} & C_{2}^{2} & C_{2}^{3} & C_{2}^{4} & C_{2}^{5} & C_{2}^{6} \\ C_{3}^{1} & C_{3}^{2} & C_{3}^{3} & C_{3}^{4} & C_{3}^{5} & C_{3}^{6} \\ C_{4}^{1} & C_{4}^{2} & C_{4}^{3} & C_{4}^{4} & C_{4}^{5} & C_{4}^{6} \end{matrix}] = [\begin{matrix} * & * & * & * & * & * \\ - \frac{11}{3} & - 3 & - \frac{5}{3} & - \frac{35}{27} & - \frac{9}{8} & - \frac{77}{75} \\ 0 & \frac{5}{9} & \frac{7}{15} & \frac{21}{50} & \frac{88}{225} & \frac{13}{35} \\ \frac{27}{28} & \frac{21}{25} & \frac{27}{35} & \frac{891}{1225} & \frac{39}{56} & \frac{33}{49} \\ \frac{80}{81} & \frac{45}{49} & \frac{55}{63} & \frac{1430}{1701} & \frac{40}{49} & \frac{1105}{1386} \end{matrix}] .

Acknowledgments

I thank Dr. John R. Rogers for encouraging me to include Section 5 and an anonymous reviewer for requesting Eqs. (6.1) and (6.2).

References and links

1. H. J. Birchall, “Lenses and their combination and arrangement in various instruments and apparatus,” U.S. patent 2,001,952 (21 May 1935).

2. H. J. Birchall, “Lens of variable focal power having surfaces of involute form,” U.S. patent 2,475,275 (7 March 1949).

3. C. W. Kanolt, “Multifocal ophthalmic lenses,” U.S. patent 2,878,721 (24 March 1959).

4. W. T. Plummer, J. G. Baker, and J. Van Tassell, “Photographic optical systems with nonrotational aspheric surfaces,” Appl. Opt. 38(16), 3572–3592 (1999). [CrossRef] [PubMed]

5. L. Wang, P. Benítez, J. C. Miñano, J. Infante, and G. Biot, “Advances in the SMS design method for imaging optics,” Proc. SPIE 8167, 81670M (2011). [CrossRef]

6. F. Muñoz, P. Benítez, and J. C. Miñano, “High-order aspherics: the SMS nonimaging design method applied to imaging optics,” Proc. SPIE 7061, 70610G, 70610G-9 (2008). [CrossRef]

7. K. H. Fuerschbach, K. P. Thompson, and J. P. Rolland, “A new generation of optical systems with φ-polynomial surfaces,” Proc. SPIE 7652, 76520C, 76520C-7 (2010). [CrossRef]

8. J. R. Rogers, “A comparison of anamorphic, keystone, and Zernike surface types for aberration correction,” Proc. SPIE 7652, 76520B, 76520B-8 (2010). [CrossRef]

9. A. Yabe, “Method to allocate freeform surfaces in axially asymmetric optical systems,” Proc. SPIE 8167, 816703, 816703-10 (2011). [CrossRef]

10. R. Steinkopf, L. Dick, T. Kopf, A. Gebhardt, S. Risse, and R. Eberhardt, “Data handling and representation of freeform surfaces,” Proc. SPIE 8169, 81690X, 81690X-9 (2011). [CrossRef]

11. P. Jester, C. Menke, and K. Urban, “B-spline representation of optical surfaces and its accuracy in a ray trace algorithm,” Appl. Opt. 50(6), 822–828 (2011). [CrossRef] [PubMed]

12. G. W. Forbes, “Manufacturability estimates for optical aspheres,” Opt. Express 19(10), 9923–9941 (2011). [CrossRef] [PubMed]

13. G. W. Forbes, “Shape specification for axially symmetric optical surfaces,” Opt. Express 15(8), 5218–5226 (2007). [CrossRef] [PubMed]

14. G. W. Forbes, “Robust, efficient computational methods for axially symmetric optical aspheres,” Opt. Express 18(19), 19700–19712 (2010). [CrossRef] [PubMed]

15. C. Zhao and J. H. Burge, “Orthonormal vector polynomials in a unit circle, Part I: Basis set derived from gradients of Zernike polynomials,” Opt. Express 15(26), 18014–18024 (2007). [CrossRef] [PubMed]

16. G. W. Forbes, “Robust and fast computation for the polynomials of optics,” Opt. Express 18(13), 13851–13862 (2010). [CrossRef] [PubMed]

Characterizing the shape of freeform optics

Abstract

1. Introduction

2. Coordinates, aperture shapes, and sag expressions

3. Orthogonalization in terms of the mean square gradient

4. Simple surface for demonstration

5. The option for a non-zero conic constant

6. Concluding remarks

Appendix A: Auxiliary polynomials

Appendix B: Additional recurrence-based processes

Acknowledgments

References and links

Cited By

Figures (8)

Equations (57)

Optics Express