1. Introduction

When either designing, fabricating, or testing optical aspheres, there are advantages to characterizing the surface’s shape in terms of orthogonal bases [1]. The coefficients in such a description can be interpreted as a spectrum-like decomposition of the shape. These individual coefficients have immediate interpretations, and they permit a part to be assessed at a glance. For {Q ^bfs _m(x)} of [1], in particular, the coefficients relate directly to the asphere’s deviation from the best-fit sphere measured along the surface normal. The fringe density in a conventional full-aperture interferogram of such a part is proportional to the rate of change of this normal departure, referred to here as the normal slope. The weighted mean square value of this normal slope is directly linked to the sum of the squares of the associated coefficients when using {Q ^bfs _m(x)}. In fact, this property served as the defining relation for the basis and facilitates design constraints when the goal is mild aspheres that are easily tested. It is important to appreciate, however, that this basis is also well suited to the characterization of arbitrarily strong aspheres. As a first step in this direction, a related design constraint was introduced for aspheres that are to be tested by using stitched interferometry [2]. More generally, a part’s manufacturability is related to the variation in the local principal curvatures of the surface, see e.g [3]. Because the normal slope couples directly to changes in the principal curvatures, this basis is well suited to a variety of manufacturability-driven constraints.

Robust and efficient computational methods are vital during the design and modeling of systems that incorporate aspheric surfaces. A trend to more complex aspheres has led to the need for more terms in the polynomial used to characterize shape. The catastrophic failure associated with round-off in these cases was recently highlighted [4] and it was shown that, when using orthogonal polynomials, recurrence relations play a vital role in avoiding this obstacle. Methods were presented in [4] for robustly computing not only the polynomials themselves, but also any linear combination of them and of their derivatives. These methods are all based on the three-term recurrence relation satisfied by standard orthogonal polynomials. It turns out, however, that {Q ^bfs _m(x)} does not satisfy a three-term recurrence. Instead, it is shown in the following sections that there is a non-standard recurrence relation that can be used to derive effective algorithms for working with this basis. In fact, the most effective algorithms are found to operate with an auxiliary set of polynomials in tandem.

2. Unconventional recurrence relations

If ρ _max denotes one half of the part’s clear aperture, a rotationally symmetric asphere’s shape is specified here in cylindrical polar coordinates as

where c is the curvature of the best-fit sphere, u is the normalized radial coordinate defined by max u:=(ρ/ρ _max), and Q_m(x) is a polynomial of order m. The coefficients written as a_m serve to characterize the asphere’s departure from its best-fit sphere. Note that I typically drop the superscript on Q ^bfs _m(x) now because this is the only one of the bases from [1] that is used in this work. This basis is constructed so that

where δ_mn is Kronecker’s delta and the normal slope of the basis members is defined by

The normalization factor in Eq. (2.2) follows from

Notice that the additive aspheric term in Eq. (2.1) corresponds to deviation measured along the axis. To first order, this is converted to departure along the surface normal by multiplying it by a cosine factor. This conversion is why the square root in the denominator in Eq. (2.1) is absent from Eq. (2.3). A sample of these basis members is plotted in Fig. 1 together with their normal slopes. The clear sine-like character of S_m(u) is intuitively consistent with the fact that the mean of the squared normal slope is just the sum of the squares of the a_m coefficients, see Eq. (15) of [1].

 

Fig. 1. A sample of the slope-orthogonal polynomials used in this work. By construction, they are not only orthogonal in slope, but normalized so that their mean square slope is unity.

Download Full Size | PDF

It is shown in Appendix A that {Q ^bfs _m(x)} is simply related to a particular family of standard polynomials. In particular, Eq. (A.12) can be written more explicitly as

where P_m(x) is a scaled Jacobi polynomial. An efficient process for evaluating the constants f_m, g_m, and h_m is given in Eqs. (A.14–16) of Appendix A. It follows from Eq. (A.4) together with 22.7.1 and 22.4.1 of [6], that these auxiliary polynomials can be generated robustly by a standard three-term recurrence relation of particularly simple form:

Equation (2.6) is initiated with P ₀(x) = 2 and P ₁(x) = 6 − 8x. A sample of these polynomials is plotted in Fig. 2. As can be seen from Eq. (A.6), the envelope drawn in Fig. 2 is given by 2u(1 − u ²), which peaks at u=3^-1/2 ≈ 0.58 where it takes the value 3^-1/2 4/3 ≈ 0.77.

Upon eliminating the P’s from Eq. (2.6) by using Eq. (2.5), an unconventional five-term recurrence relation emerges for {Q ^bfs _m(x)}. It turns out that, rather than doing this explicitly, it is more effective for several reasons to operate with {P_m(x)} and {Q ^bfs _m(x)} in tandem. First, notice that {Q ^bfs _m(x)} can be generated robustly to arbitrary orders upon rewriting Eq. (2.5) as an unconventional three-term recurrence relation:

with Q ₀(x) = 1 and Q ₁(x) = 19^−1/2 (13 − 16x). More importantly, as shown in the next section, the change of basis matrix determined in Appendix A can be used to great effect for a number of critical operations when working with aspheres in these terms.

 

Fig. 2. A sample of the scaled Jacobi polynomials used in this work. The plot at left is for comparison with Fig. 1, and the one on the right holds up to m = 20 to reveal their envelope.

Download Full Size | PDF

3. Linear combinations and their derivatives

The key piece of interest in Eq. (2.1) can n be written as S(u ²), where S is defined by

If the list of these basis elements is denoted by q(x), then S(x) = a·q(x). Similarly, write the auxiliary basis, i.e. {P_m(x)}, as p(x), and define a second function, say T(x) := b·p(x). Given that p(x) = Lq(x) where L is the lower-triangular band matrix of Eq. (2.5) or Eq. (A.12), it follows that T(x) = b·[Lq(x)] = (L^T b)·q(x). It is now evident that S(x) ≡ T(x) when a = L^T b. As a consequence, we can change basis from {P_m(x)} to {Q ^bfs _m(x)} simply by using

for m = 0, 1, 2,… M − 2 and then

The inverse of this transformation is also easily realized by reversing these steps: start with

and, for m = M − 2 down to 0, use

The recurrence relation in Eq. (3.5) is clearly just a re-arranged version of Eq. (3.2). In this way, with about 5M arithmetic operations, it is possible to interchange between these two bases. [Such a process typically involves O(M ²) operations, but the band matrix accelerates it in this case.]

An advantage of swapping bases is that, as shown in [4], Clenshaw’s method [7] allows

to be computed robustly by exploiting the recurrence relation in Eq. (2.6): Simply start with

and work down with

from m = M − 2 to 0. The desired end result is then given —see Eq. (A.9) of [4] — as

The final step follows from the expressions for the lowest two basis members that were given after Eq. (2.6). That is, with the roughly 5M arithmetic operations needed to evaluate the α’s via Eqs. (3.7) and (3.8), S(x) can be computed without evaluating a single member of either basis.

It follows from Eqs. (3.7) and (3.8) that α_m is a polynomial of order M − m in x. As discussed in [4], Smith [8] demonstrated that any derivative of S(x) of Eq. (3.6) can be evaluated with similar ease. If the j’th derivative is written as S ^(j)(x), then start with

and work down from m = M − 2 to 0 by using

The desired end result is then

Note that this process assumes that the (j − 1)’th derivative has been evaluated before going after the j’th, and that α ⁽⁰⁾ _m is equal to α_m. It follows from Eq. (2.1) that, with ϕ:=(1 − c ² ρ ²)^1/2, the surface’s sag and its first two derivatives are given in these terms by

In typical applications, the sag of a specific part and its derivatives will be evaluated for multiple values of ρ. In such cases, {a_m} is converted to {b_m} only once, and these serve as the input —along with various ρ values— to the methods just described. For example, to plot a specific member of {Q ^bfs _m(x)} or its derivatives, rather than use a process like that discussed at Eq. (2.7), it is better to set all of {a_m} to zero except the one of interest and convert this list to the associated {b_m} via Eqs. (3.4) and (3.5). The desired results then follow efficiently from Eqs. (3.7) and (3.9) and (3.10) to (3.12). It is also worth noting that the part’s axial curvature follows upon taking ρ = 0 (hence u = 0 and ϕ = 1) in Eq. (3.15):

This final expression follows from the fact that P_m(0) = 2(2m + 1), which can be derived inductively from Eq. (2.6). By either solving Eq. (3.16) for c, or leaving it in terms of c _axial, it is possible to consider the independent curvature variable to be either the axial curvature or the best-fit curvature when using this representation. Taking c _axial to be the more fundamental entity can be helpful for paraxial analysis and constraints in the design process.

4. Fitting to this orthogonal basis

Consider an aspheric surface specified by a general sag function of the form

where f(ρ) is symmetric and f(0) = 0. To express this in the form given in Eq. (2.1), the first step is to determine the curvature of the best-fit sphere by using the requirement that it meets the asphere at both ρ = 0 and ρ = ρ _max :

The task is then reduced to choosing {a_m} so that

Since {a_m} enters this relation linearly, a standard linear least-squares process can determine the solution via an (M + 1) × (M + 1) matrix inversion. This is adequate for most purposes. It is impressive, however, that swapping to the auxiliary basis —where a and Q in Eq. (4.3) are replaced by b and P — opens an option for an elegant and more efficient solution.

After re-arranging Eq. (4.3) and inserting ρ = uρ _max, the task in terms of the auxiliary basis is to determine {b_m} so that

The last relation serves to define F(u) which, provided f′(0) = 0, remains well defined even at the endpoints of 0 ≤ u ≤ 1: the term inside the braces has zeros that remove the singularities that are potentially caused by the denominator outside the braces. It follows from Eq. (B.7) of Appendix B that a particular solution for these coefficients is just

where

As discussed in Appendix B, Eq. (4.5) is a standard discrete cosine transform (namely DCT-IV) and, when more than just a few terms are required in the fit, it may be better implemented via FFT-like options. In our context, it is unusual for the spectrum to spill significantly beyond just tens of terms. For most purposes, therefore, troubles with aliasing are insignificant if we use Eq. (4.5) with N about 16 or 32. Remember that, as shown in Fig. 2, the maximum contribution to the normal departure from the best-fit sphere associated with any b_m is about 0.8b_m. It is therefore straightforward to truncate the list of these coefficients so that the magnitude of the dropped terms is less than whatever cut-off is desired. As a final step, Eqs. (3.2) and (3.3) yield {a_m}.

For the purpose of demonstration, consider a parabola specified by c _axial = 20mm and ρ _max = 20mm. It follows from Eq. (4.2) that the best-fit curvature is c = 25mm. A scaled lens drawing and a plot of the associated sag-based departure from best-fit sphere are presented in Fig. 3. The fit coefficients that result from Eq. (4.5) with N = 32 for this case are plotted in Fig. 4. As it does here, the magnitude of the coefficients generally decays rapidly with order. Notice that these coefficients hit the floor of numerical round-off for double precision at about m = 20. The order at which this happens depends on the characteristics of the shape. For nm-level accuracy, since the maximum contribution from each term is about 0.8 b_m, this list of coefficients can be truncated at about the level of the solid gray line in the plot.

 

Fig. 3. A cross-section of the parabola used for demonstration (in blue) is shown with its bestfit sphere (dotted red). The aspheric departure (viz. sag difference) is plotted at right.

Download Full Size | PDF

To facilitate code verification, the eight coefficients above the solid gray line in Fig. 4 are found to be

If the last term in this subset is also dropped, the fit error can therefore be expected to resemble the m = 7 plot of Fig. 2 with a scale factor of 2.6nm. This is perfectly consistent with the computed fit error plotted at right in Fig. 4. After converting the seven retained coefficients to {a_m} according to Eqs. (3.2) and (3.3), the specification for this surface is found to be

It was these rounded coefficients that were used to compute the fit error plotted in Fig. 4. Interestingly, for any N ≥ 16, just the last few digits vary in only the smaller members of the list in Eq. (4.7). Even with N = 8 in this case, the aliasing effects are so small that nm-level accuracy is achieved in the resulting fit.

 

Fig. 4. Log plot of {b_m} in nm units. The dots are red/blue for positive/negative values. The 1nm reference level is drawn as a solid gray line. The plot at right is the error in the fit that results if all but the last red dot above the 1nm line are retained (i.e. m = 0 to 6 are kept).

Download Full Size | PDF

5. Annular apertures

For obstructed mirror systems, the aspheric departure of each surface is defined (and measured) over only the appropriate annulus, say ερ _max < ρ < ρ _max where 0 < ε < 1. With a minor compromise, the ideas described above can be applied to such cases with minimal change. The best-fit sphere is now chosen to pass through both the inner and the outer edge of the annulus, and Eq. (2.1) is replaced by

This turns out to be orthonormal in slope if we choose the inner product for the mean square slope in this case to be

where Eq. (5.1) means that the normal slope is now defined by

The normalization factor in Eq. (5.2) follows from

Orthonormality can be confirmed by changing variables to x:=(u ² − ε ²)/(1 − ε ²) in order to see that Eqs. (5.2) and (5.3) reduce to precisely the inner product in Eq. (A.5) of Appendix A.

 

Fig. 5. The first ten basis members and their slopes for an obstructed aperture with ε = 0.25.

Download Full Size | PDF

The compromise here is that the weight function in Eq. (5.2) goes to zero at the inner edge of the annulus, i.e. at u = ε. As a result, as is evident in Fig. 5, the slope maps are not as sine-like as those in Fig. 1, but tend to grow at this edge. Ideally, a new set of polynomials would be defined by replacing the weight in Eq. (5.2) by a weight that diverges similarly at each edge, namely u[(u ² − ε ²)(1 − u ²)]^−1/2. In my opinion, this complication is unjustified, and it is better to re-use the results derived above for the unobstructed case with minor adaptation. The key is that, at its heart, Eq. (5.1) is built around a sum like that in Eq. (3.1), hence also Eq. (3.6). On account of the factor 2/x ^1/2 in Eq. (A.6), the envelope in the analogue of Fig. 2 is now 2(1 − u ²)[(u ² − ε ²)/(1 − ε)]^1/2/(1 + ε), with a peak value of $\frac{4}{3}\left(1-\epsilon \right){\left[\frac{\left(1+\epsilon \right)}{3}\right]}^{\frac{1}{2}}$ . Of course, the DCT is now applied to the resulting analog of Eq. (4.4). In the limit of small ε, all these results reduce precisely to those for the unobstructed case.

6. Concluding remarks

Although two bases are used in this work, the end user need only ever consider one of them; the other can be hidden within the software as an internal computational aid. Because manufacturability is more closely coupled to the normal slope than to sag, I regard {Q ^bfs _m(x)} as the principal basis, i.e. I prefer to think in terms of a rather than b. If the illuminated region plus whatever margin is required (for fabrication) differs significantly from the disc of radius ρ _max, however, the constraints based on a will lose validity. It can therefore be helpful during optimization to be able to adjust an asphere’s clear aperture while preserving its shape. Fortunately, the method described in Section 4 allows this adjustment to be performed efficiently. The elegance and power of all the methods described above mean that such processes may well spread into the domain of characterizing mid-spatial frequencies on optical surfaces: with Eqs. (4.4)–(4.6), it is just as easy now to use hundreds of terms to fit a measured part. High orders can also be valuable in simulations during design for tolerancing, etc. Although it is beyond our current scope, it is natural for these sorts of applications —as well as for the growing field of freeform optics— to include non-rotationally symmetric terms by using a Fourier-series-based process that is similar to that used in the construction of the Zernike polynomials.

The methods reported in this work lead to concise code that offers a remarkably efficient manner to work with aspheres in terms of {Q ^bfs _m(x)}. These methods are robust to arbitrary orders and facilitate characterizations like Eq. (4.8) that have been found to hold typically three times fewer digits than the traditional representation [11]. Their spectrum-like properties also allow the manufacturing challenge to be estimated at a glance. For the example in Sec. 4, for instance, the rapid decay in the coefficients reveals that the shape is dominated by the lowest order term, i.e. the red curve of Fig. 1. What’s more, its peak value is therefore roughly 0.25 times the first coefficient of the list in Eq. (4.8). Aspheric departure plots —in this case, the solid red curve in Fig. 3— can typically be anticipated in this way with no computation. What’s more, the number of terms that are essential becomes immediately self-evident. Perhaps most important of all, these new processes avoid the catastrophic round-off failure that has plagued this field. Whether in the context of design, fabrication, or testing, the benefits of an orthogonal basis are compelling, and these new recurrence-based methods mean that any additional operational costs are insignificant.