## Robust, efficient computational methods for axially symmetric optical aspheres |

Optics Express, Vol. 18, Issue 19, pp. 19700-19712 (2010)

http://dx.doi.org/10.1364/OE.18.019700

Acrobat PDF (898 KB)

### Abstract

Whether in design or the various stages of fabrication and testing, an effective representation of an asphere’s shape is critical. Some algorithms are given for implementing tailored polynomials that are ideally suited to these needs. With minimal coding, these results allow a recently introduced orthogonal polynomial basis to be employed to arbitrary orders. Interestingly, these robust and efficient methods are enabled by the introduction of an auxiliary polynomial basis.

© 2010 Optical Society of America

## 1. Introduction

4. G. W. Forbes, “Robust and fast computation for the polynomials of optics,” Opt. Express **18**, 13851–13862, (2010) http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-18-13-13851. [CrossRef] [PubMed]

4. G. W. Forbes, “Robust and fast computation for the polynomials of optics,” Opt. Express **18**, 13851–13862, (2010) http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-18-13-13851. [CrossRef] [PubMed]

*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

*ρ*

_{max}denotes one half of the part’s clear aperture, a rotationally symmetric asphere’s shape is specified here in cylindrical polar coordinates as

*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*serve to characterize the asphere’s departure from its best-fit sphere. Note that I typically drop the superscript on

_{m}*Q*

^{bfs}

_{m}(

*x*) now because this is the only one of the bases from [1

1. G. W. Forbes, “Shape specification for axially symmetric optical surfaces,” Opt. Express **15**, 5218–5226, (2007) http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-15-8-5218. [CrossRef] [PubMed]

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

_{mn}*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*coefficients, see Eq. (15) of [1

_{m}1. G. W. Forbes, “Shape specification for axially symmetric optical surfaces,” Opt. Express **15**, 5218–5226, (2007) http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-15-8-5218. [CrossRef] [PubMed]

*P*(

_{m}*x*) is a scaled Jacobi polynomial. An efficient process for evaluating the constants

*f*,

_{m}*g*, and

_{m}*h*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:

_{m}*P*

_{0}(

*x*) = 2 and

*P*

_{1}(

*x*) = 6 − 8

*x*. 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 2

*u*(1 −

*u*

^{2}), which peaks at

*u*=3

^{-1/2}≈ 0.58 where it takes the value 3

^{-1/2}4/3 ≈ 0.77.

*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:

*Q*

_{0}(

*x*) = 1 and

*Q*

_{1}(

*x*) = 19

^{−1/2}(13 − 16

*x*). 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.

## 3. Linear combinations and their derivatives

**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*) = L

**q**(

*x*) where L is the lower-triangular band matrix of Eq. (2.5) or Eq. (A.12), it follows that

*T*(

*x*) =

**b**·[L

**q**(

*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

*m*= 0, 1, 2,…

*M*− 2 and then

*m*=

*M*− 2 down to 0, use

*M*arithmetic operations, it is possible to interchange between these two bases. [Such a process typically involves

*O*(

*M*

^{2}) operations, but the band matrix accelerates it in this case.]

4. G. W. Forbes, “Robust and fast computation for the polynomials of optics,” Opt. Express **18**, 13851–13862, (2010) http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-18-13-13851. [CrossRef] [PubMed]

7. C. W. Clenshaw, “A Note on the Summation of Chebyshev Series,” Math. Tables Other Aids Comput. **9**, 118–120 (1955), http://www.jstor.org/stable/2002068.

*m*=

*M*− 2 to 0. The desired end result is then given —see Eq. (A.9) of [4

**18**, 13851–13862, (2010) http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-18-13-13851. [CrossRef] [PubMed]

*M*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.

*α*is a polynomial of order

_{m}*M*−

*m*in

*x*. As discussed in [4

**18**, 13851–13862, (2010) http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-18-13-13851. [CrossRef] [PubMed]

8. F. J. Smith, “An Algorithm for Summing Orthogonal Polynomial Series and their Derivatives with Applications to Curve-Fitting and Interpolation,” Math. Comput. **19**, 33–36 (1965). [CrossRef]

*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

*m*=

*M*− 2 to 0 by using

*j*− 1)’th derivative has been evaluated before going after the

*j*’th, and that

*α*

^{(0)}

_{m}is equal to

*α*. It follows from Eq. (2.1) that, with

_{m}*ϕ*:=(1 −

*c*

^{2}

*ρ*

^{2})

^{1/2}, the surface’s sag and its first two derivatives are given in these terms by

*ρ*. In such cases, {

*a*} is converted to {

_{m}*b*} only once, and these serve as the input —along with various

_{m}*ρ*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*} to zero except the one of interest and convert this list to the associated {

_{m}*b*} 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

_{m}*ρ*= 0 (hence

*u*= 0 and

*ϕ*= 1) in Eq. (3.15):

*P*(0) = 2(2

_{m}*m*+ 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

*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}:

*a*} so that

_{m}*a*} enters this relation linearly, a standard linear least-squares process can determine the solution via an (

_{m}*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.

*ρ*=

*uρ*

_{max}, the task in terms of the auxiliary basis is to determine {

*b*} so that

_{m}*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

*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*is about 0.8

_{m}*b*. 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 {

_{m}*a*}.

_{m}*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*, this list of coefficients can be truncated at about the level of the solid gray line in the plot.

_{m}*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*} according to Eqs. (3.2) and (3.3), the specification for this surface is found to be

_{m}*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.

## 5. Annular apertures

*ερ*

_{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

*x*:=(

*u*

^{2}−

*ε*

^{2})/(1 −

*ε*

^{2}) in order to see that Eqs. (5.2) and (5.3) reduce to precisely the inner product in Eq. (A.5) of Appendix A.

*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*

^{2}−

*ε*

^{2})(1 −

*u*

^{2})]

^{−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*

^{2})[(

*u*

^{2}−

*ε*

^{2})/(1 −

*ε*)]

^{1/2}/(1 +

*ε*), with a peak value of

*ε*, all these results reduce precisely to those for the unobstructed case.

## 6. Concluding remarks

*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

11. G. W. Forbes, “Can you make/measure this asphere for me?”, Frontiers in Optics, OSA Technical Digest (2009), http://www.opticsinfobase.org/abstract.cfm?URI=FiO-2009-FThH1.

## Appendix A: Recurrence involving a Jacobi polynomial

*u*to

*x*:=

*u*

^{2}and using Eq. (2.3) to see that this condition is equivalent to

*ϕ*(

_{n}*x*) is any polynomial of order

*n*, the Jacobi polynomials satisfy [5

5. E. W. Weisstein, “Jacobi Polynomial” from MathWorld—A Wolfram Web Resource. http://mathworld.wolfram.com/JacobiPolynomial.html, see esp. Equations (10–14).

5. E. W. Weisstein, “Jacobi Polynomial” from MathWorld—A Wolfram Web Resource. http://mathworld.wolfram.com/JacobiPolynomial.html, see esp. Equations (10–14).

*Q*(

_{m}*x*) is replaced by

*P*

^{(−1/2,1/2)}

_{m}(2

*x*− 1) throughout Eq. (A.1), the expression on its left-hand side vanishes whenever ∣

*m*−

*n*∣ > 2. That is, once expanded out, the integral of each term is found to be zero, so the associated matrix is quindiagonal: it has the main diagonal and two bands of non-zero elements on either side. This reveals a useful link between this particular family of Jacobi polynomials and

*Q*

^{bfs}

_{m}(

*x*). It turns out to be convenient to express what follows in terms of

*n*!! = (

*n*− 2)!!

*n*, and 0!! = (−1)!! = 1, so 7!! = 7·5·3·1 etc.

*Q*> =

_{m}, Q_{n}*δ*. More interestingly, it was established in the previous paragraph that <

_{mn}*P*> is quindiagonal. This can be seen more directly upon using an explicit expression for

_{m}, P_{n}*P*(

_{m}*x*) that follows from a link to the Chebyshev polynomials of the first kind, see 22.3.15 and 22.5.29 of [6], namely

*x*= cos

^{2}

*θ*, it follows from Eqs. (A.5) and (A.6) that

*P*,

_{m}*P*> is quindiagonal. What’s more, it becomes straightforward to evaluate its elements. The three independent non-zero bands in this symmetric matrix are next handled separately.

_{n}*H*:= <

_{m}*P*

_{m+2},

*P*>, it follows that the only non-zero contribution in Eq. (A.7) is of the form

_{m}*G*:= <

_{m}*P*

_{m+1},

*P*>, it follows similarly (although now with two non-zero terms) that

_{m}*F*:= <

_{m}*P*,

_{m}*P*>, it is found that the three non-zero terms give

_{m}*m*= 0, however, the cos[(2

*m*+ 1)

*θ*] and cos[(2

*m*− 1)

*θ*] terms are no longer orthogonal. Instead, they add directly and it is readily seen that

*F*

_{0}= 4. As shown next, these matrix elements play a fundamental role in the relation between {

*Q*

^{bfs}

_{m}(

*x*)} and {

*P*(

_{m}*x*)}.

*n*’th order polynomial can obviously be expressed as a combination of

*Q*(

_{m}*x*) for

*m*= 0,1,…

*n*means that this matrix is lower-triangular, i.e.

*L*= 0 for

_{ij}*j*>

*i*. By using <

*Q*,

_{m}*Q*> =

_{n}*δ*, it now follows that

_{mn}*P*,

_{m}*P*> is just the product of the matrix with elements

_{n}*L*and its transpose, say LL

_{jk}^{T}. It follows that L is also a band matrix: the only non-zero elements are along the diagonal and the two bands below the diagonal.

*f*:=

_{m}*L*,

_{mm}*g*:=

_{m}*L*

_{m+1m}, and

*h*:=

_{m}*L*

_{m+2m}, these elements can be found by starting with

*f*

_{0}= 2,

*f*

_{1}= 19

^{1/2}/2,

*g*

_{0}=−1/2 and working up from

*m*= 2 by applying

*m*’th row of this matrix from the earlier rows, and this process is robust to arbitrary orders. It turns out that, for large

*m*,

*f*≈ −

_{m}*h*≈ 2

_{m}^{−1/2}

*m*and

*g*≈ −2

_{m}^{−1/2}.

## Appendix B: Determining a fit in terms of the auxiliary Jacobi polynomial

*g*(

*x*), it follows that the fit that minimizes the mean square error defined by

*b*can be found by averaging adjacent Fourier coefficients from the Fourier series of the periodic and even function

_{m}*g*[(1 + cos

*ψ*)/2]. Further, because cosines are orthogonal over summation as well as integration, this opens striking new options. In particular, when only a finite number, say

*N*, of these Fourier coefficients are non-zero, Eq. (B.5) can be computed exactly by using the discrete cosine transform (DCT) involving a sum over

*N*points. Of course, the DCT can be accelerated by using the FFT in a variety of ways, see Sec 12.3 of [9].

*π*/(2

*N*) in one of two obvious ways:

10. A useful overview is given at http://en.wikipedia.org/wiki/Discrete_cosine_transform.

*N*.

## References and links

1. | G. W. Forbes, “Shape specification for axially symmetric optical surfaces,” Opt. Express |

2. | G. W. Forbes and C. P. Brophy, “Designing cost-effective systems that incorporate high-precision aspheric optics”, SPIE Optifab (2009) TD06–25 (1). Available at http://www.qedmrf.com. |

3. | C. du Jeu, “Criterion to appreciate difficulties of aspherical polishing,” Proc. SPIE |

4. | G. W. Forbes, “Robust and fast computation for the polynomials of optics,” Opt. Express |

5. | E. W. Weisstein, “Jacobi Polynomial” from MathWorld—A Wolfram Web Resource. http://mathworld.wolfram.com/JacobiPolynomial.html, see esp. Equations (10–14). |

6. | M. Abramowitz and I. Stegun, |

7. | C. W. Clenshaw, “A Note on the Summation of Chebyshev Series,” Math. Tables Other Aids Comput. |

8. | F. J. Smith, “An Algorithm for Summing Orthogonal Polynomial Series and their Derivatives with Applications to Curve-Fitting and Interpolation,” Math. Comput. |

9. | W. Press, S. Teukolsky, W. Vetterling, and B. Flannery, |

10. | A useful overview is given at http://en.wikipedia.org/wiki/Discrete_cosine_transform. |

11. | G. W. Forbes, “Can you make/measure this asphere for me?”, Frontiers in Optics, OSA Technical Digest (2009), http://www.opticsinfobase.org/abstract.cfm?URI=FiO-2009-FThH1. |

**OCIS Codes**

(000.4430) General : Numerical approximation and analysis

(220.1250) Optical design and fabrication : Aspherics

(220.4610) Optical design and fabrication : Optical fabrication

(220.4830) Optical design and fabrication : Systems design

(220.4840) Optical design and fabrication : Testing

**ToC Category:**

Optical Design and Fabrication

**History**

Original Manuscript: June 1, 2010

Revised Manuscript: July 12, 2010

Manuscript Accepted: July 12, 2010

Published: September 1, 2010

**Citation**

G. W. Forbes, "Robust, efficient computational methods for axially symmetric optical aspheres," Opt. Express **18**, 19700-19712 (2010)

http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-18-19-19700

Sort: Year | Journal | Reset

### References

- G. W. Forbes, "Shape specification for axially symmetric optical surfaces," Opt. Express 15, 5218-5226 (2007), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-15-8-5218. [CrossRef] [PubMed]
- G. W. Forbes, and C. P. Brophy, "Designing cost-effective systems that incorporate high-precision aspheric optics," SPIE Optifab (2009) TD06-25 (1). Available at http://www.qedmrf.com.
- C. du Jeu, "Criterion to appreciate difficulties of aspherical polishing," Proc. SPIE 5494, 113-121 (2004), doi:10.1117/12.551420. [CrossRef]
- G. W. Forbes, "Robust and fast computation for the polynomials of optics," Opt. Express 18, 13851-13862 (2010), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-18-13-13851. [CrossRef] [PubMed]
- E. W. Weisstein, "Jacobi Polynomial" from MathWorld—A Wolfram Web Resource. http://mathworld.wolfram.com/JacobiPolynomial.html, see esp. Eqs. (10)-(14).
- M. Abramowitz, and I. Stegun, Handbook of Mathematical Functions (Dover, 1978), Chap. 22.
- C. W. Clenshaw, "A Note on the Summation of Chebyshev Series," Math. Tables Other Aids Comput. 9, 118-120 (1955), http://www.jstor.org/stable/2002068.
- F. J. Smith, "An Algorithm for Summing Orthogonal Polynomial Series and their Derivatives with Applications to Curve-Fitting and Interpolation," Math. Comput. 19, 33-36 (1965). [CrossRef]
- W. Press, S. Teukolsky, W. Vetterling, and B. Flannery, Numerical Recipes: The Art of Scientific Computing (Cambridge University Press, 1992) Section 2.9.
- A useful overview is given at http://en.wikipedia.org/wiki/Discrete_cosine_transform.
- G. W. Forbes, "Can you make/measure this asphere for me?", Frontiers in Optics, OSA Technical Digest (2009), http://www.opticsinfobase.org/abstract.cfm?URI=FiO-2009-FThH1.

## Cited By |
Alert me when this paper is cited |

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

« Previous Article | Next Article »

OSA is a member of CrossRef.