OSA's Digital Library

Optics Express

Optics Express

  • Editor: C. Martijn de Sterke
  • Vol. 18, Iss. 19 — Sep. 13, 2010
  • pp: 19700–19712
« Show journal navigation

Robust, efficient computational methods for axially symmetric optical aspheres

G.W. Forbes  »View Author Affiliations


Optics Express, Vol. 18, Issue 19, pp. 19700-19712 (2010)
http://dx.doi.org/10.1364/OE.18.019700


View Full Text Article

Acrobat PDF (898 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

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

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

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]

] and it was shown that, when using orthogonal polynomials, recurrence relations play a vital role in avoiding this obstacle. Methods were presented in [4

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]

] 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

z=cρ21+1c2ρ2+u2(1u2)1c2ρ2m=0MamQm(u2),
(2.1)

where c is the curvature of the best-fit sphere, u is the normalized radial coordinate defined by max u:=(ρ/ρ max), and Qm(x) is a polynomial of order m. The coefficients written as am 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

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]

] that is used in this work. This basis is constructed so that

(2π)01Sm(u)Sn(u)11u2du=δmn,
(2.2)

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

Sm(u)ddu[u2(1u2)Qm(u2)].
(2.3)

The normalization factor in Eq. (2.2) follows from

0111u2du=π2.
(2.4)

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 Sm(u) is intuitively consistent with the fact that the mean of the squared normal slope is just the sum of the squares of the am coefficients, see Eq. (15) of [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]

].

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.

Pm(x)=fmQm(x)+gm1Qm1(x)+hm2Qm2(x),
(2.5)

where Pm(x) is a scaled Jacobi polynomial. An efficient process for evaluating the constants fm, gm, and hm 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

6. M. Abramowitz and I. Stegun, Handbook of Mathematical Functions (Dover, 1978), Chap. 22.

], that these auxiliary polynomials can be generated robustly by a standard three-term recurrence relation of particularly simple form:

Pm+1(x)=(24x)Pm(x)Pm1(x).
(2.6)

Equation (2.6) is initiated with P 0(x) = 2 and P 1(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 2), 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 {Pm(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:

Qm+1(x)=[Pm+1(x)gmQm(x)hm1Qm1(x)]fm+1,
(2.7)

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.

3. Linear combinations and their derivatives

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

S(x)m=0MamQm(x).
(3.1)

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. {Pm(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)] = (LT bq(x). It is now evident that S(x) ≡ T(x) when a = LT b. As a consequence, we can change basis from {Pm(x)} to {Q bfs m(x)} simply by using

am=fmbm+gmbm+1+hmbm+2,
(3.2)

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

aM1=fM1bM1+gM1bMandaM=fMbM.
(3.3a,b)

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

bM=aMfMandbM1=(aM1gM1bM)fM1,
(3.4a,b)

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

bm=(amgmbm+1hmbm+2)fm.
(3.5)

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 2) operations, but the band matrix accelerates it in this case.]

An advantage of swapping bases is that, as shown in [4

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]

], Clenshaw’s method [7

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.

] allows

S(x)m=0MbmPm(x),
(3.6)

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

αM=bMandαM1=bM1+(24x)αM,
(3.7a,b)

and work down with

αm=bm+(24x)αm+1αm+2,
(3.8)

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

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]

] — as

S(x)=[α0(24x)α1]P0(x)+α1P1(x)
=2(α0+α1).
(3.9)

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 Mm in x. As discussed in [4

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]

], Smith [8

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]

] 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

αMj+1(j)=0andαMj(j)=4jαMj+1(j1),
(3.10a,b)

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

αm(j)=(24x)αm+1(j)αm+2(j)4jαm+1(j1).
(3.11)

The desired end result is then

S(j)(x)=2(α0(j)+α1(j)).
(3.12)

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

z(ρ)=cρ21+φ+u2(1u2)φS(u2),
(3.13)
z(ρ)=cρφ+u[1+φ2u2(1+3φ2)]ρmaxφ3S(u2)+2u3(1u2)ρmaxφS(u2),
(3.14)
z(ρ)=cφ3+3φ23u2(1+φ2+2φ4)ρmax2φ5S(u2)
+2u2[2+3φ2u2(2+7φ2)]ρmax2φ3S(u2)+4u4(1u2)ρmax2φS(u2).
(3.15)

In typical applications, the sag of a specific part and its derivatives will be evaluated for multiple values of ρ. In such cases, {am} is converted to {bm} 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 {am} to zero except the one of interest and convert this list to the associated {bm} 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):

caxial=z(0)=c+2ρmax2S(0)=c+4ρmax2m=0M(2m+1)bm.
(3.16)

This final expression follows from the fact that Pm(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

z=f(ρ),
(4.1)

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 :

c=2f(ρmax)ρmax2+f(ρmax)2.
(4.2)

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

f(ρ)cρ21+1c2ρ2+u2(1u2)1c2ρ2m=0MamQm(u2).
(4.3)

Since {am} 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 ρ = max, the task in terms of the auxiliary basis is to determine {bm} so that

m=0MbmPm(u2)1(ucρmax)2u2(1u2){f(uρmax)c(uρmax)21+1(ucρmax)2}F(u).
(4.4)

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

bm(1)m2Nj=0N1Fjcos[πN(m+12)(j+12)],
(4.5)

where

Fjcos(π2N(j+12))F[cos(π2N(j+12))].
(4.6)

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 bm is about 0.8bm. 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 {am}.

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 bm, 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.

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

b={1009010.04959,2770.64974485,4739.30847163
1172.09704743,257.270488293,55.4172061289,
11.966650385,2.60463667585}nm.
(4.7)

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

a={2019004,7143,13944,4190,1095,283,68}nm.
(4.8)

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 {bm} 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).

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

z=cρ21+1c2ρ2+(u2ε2)(1u2)(1ε2)(1+ε)(1c2ρ2)m=0MamQm(u2ε21ε2),
(5.1)

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

[2(1ε)π]ε1Sm(u)Sn(u)1uu2ε21u2du,
(5.2)

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

Sm(u)1(1ε2)1+εddu[(u2ε2)(1u2)Qm(u2ε21ε2)].
(5.3)

The normalization factor in Eq. (5.2) follows from

ε11uu2ε21u2du=(1ε)π2.
(5.4)

Orthonormality can be confirmed by changing variables to 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.

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

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 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 43(1ε)[(1+ε)3]12 . 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

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.

]. 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.

Appendix A: Recurrence involving a Jacobi polynomial

The integral in Eq. (2.2) can be re-written by changing the variable of integration from u to x := u 2 and using Eq. (2.3) to see that this condition is equivalent to

(2π)01{[(12x)Qm(x)+x(1x)Qm(x)]
×[(12x)Qn(x)+x(1x)Qn(x)]}4x1xdx=δmn,
(A.1)

where primes denote derivatives. If ϕ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).

]

01φn(x)Pm(α,β)(2x1)(1x)αxβdx=0,form>n,
(A.2)

and

ddx[Pn(α,β)(x)]=12(n+α+β+1)Pn1(α+1,β+1)(x).
(A.3)

Pm(x)(1)m2(2m)!!(2m1)!!Pm(12,12)(2x1),
(A.4)

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

As suggested by Eqs. (2.2) and (2.3) and the change of variables used for Eq. (A.1), the natural inner-product for this context can be defined as

f,g(2π)01ddx[x(1x)f(x)]ddx[x(1x)g(x)]4x1xdx.
(A.5)

Pm(x)=(1)m2cos[2m+12arccos(2x1)]x,
(A.6)

Upon changing the variable of integration by using x = cos2 θ, it follows from Eqs. (A.5) and (A.6) that

Pm,Pn=(2π)0π2Tm(θ)Tn(θ)dθ=(1π)π2π2Tm(θ)Tn(θ)dθ,
(A.7)

where

Tm(θ)(1)m2sinθddθ{cosθsin2θcos[(2m+1)θ]}
=(1)m{(m+2)cos[(2m+3)θ]+cos[(2m+1)θ](m1)cos[(2m1)θ]}.
(A.8)

First, with Hm := <P m+2, Pm>, it follows that the only non-zero contribution in Eq. (A.7) is of the form

Hm=(1π)(m+2)(m+1)π2π2cos2[(2m+3)θ]dθ
=(m+2)(m+1)2.
(A.9)

Next, for Gm := <P m+1, Pm>, it follows similarly (although now with two non-zero terms) that

Gm=(1π){(m+2)π2π2cos2[(2m+3)θ]dθmπ2π2cos2[(2m+1)θ]dθ}
=1,
(A.10)

and finally, for Fm := <Pm, Pm>, it is found that the three non-zero terms give

Fm=(1π)[(m+2)2+1+(m1)2]π2
=(m2+m+3),form>0.
(A.11)

When m = 0, however, the cos[(2m + 1)θ] and cos[(2m − 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 {Pm(x)}.

Pi(x)=jLijQj(x).
(A.12)

The fact that any n’th order polynomial can obviously be expressed as a combination of Qm(x) for m = 0,1,… n means that this matrix is lower-triangular, i.e. Lij = 0 for j > i. By using <Qm, Qn> = δmn, it now follows that

Pm,Pn=jLmjQj(x),kLnkQk(x)
=jkLmjLnkQj(x),Qk(x)=jLmjLnj.
(A.13)

That is, the quindiagonal, <Pm, Pn> is just the product of the matrix with elements Ljk and its transpose, say LLT. 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.

The elements of the change-of-basis matrix, i.e. L, can now be evaluated by using the standard process of Cholesky decomposition [9

9. W. Press, S. Teukolsky, W. Vetterling, and B. Flannery, Numerical Recipes: The Art of Scientific Computing (Cambridge University Press, 1992) Section 2.9.

]. That is, if the non-zero elements of the change-of-basis matrix are denoted by fm := Lmm, gm := L m+1m, and hm := L m+2m, these elements can be found by starting with f 0 = 2, f 1 = 191/2/2, g 0 =−1/2 and working up from m = 2 by applying

hm2=Hm2fm2=m(m1)(2fm2),
(A.14)
gm1=(Gm1gm2hm2)fm1=(1+gm2hm2)fm1,
(A.15)
fm=Fmgm12hm22=m(m+1)+3gm12hm22,
(A.16)

in this order. The first equality in each of Eqs. (A14–16) is just for explanatory purposes; the second expression is all that is required. For reference, the initial sub-block of this change-of-basis matrix is

[f000000g0f10000h0g1f20000h1g2f30000h2g3f40000h3g4f5]=[20000012194000012521941019000061917219012509100000321910912509062595090000201050947321318311225607259].
(A.17)

Equations (A.14) and (A.16) yield the 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, fm ≈ −hm ≈ 2−1/2 m and gm ≈ −2−1/2.

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

The auxiliary polynomials used here are defined by Eq. (A.4) and satisfy

12π01Pm(x)Pm(x)x1xdx=δmn.
(B.1)

For any given g(x), it follows that the fit that minimizes the mean square error defined by

E212π01{g(x)mbmPm(x)}2x1xdx,
(B.2)

is given simply by

bm=12π01g(x)Pm(x)x1xdx.
(B.3)

Upon changing the variable of integration to x = cos2 θ, Eqs. (A.6) and (B.3) lead to

bm=(1)mππ2π2g(cos2θ)cosθcos[(2m+1)θ]dθ.
(B.4)

That is, each of these coefficients is given uniquely and explicitly by a simple integral.

There are various options for efficiently evaluating the integral in Eq. (B.4). One follows upon using elementary trigonometric identities and rescaling the variable of integration:

bm=(1)m2ππ2π2g(cos2θ){cos[2(m+1)θ]+cos[2mθ]}dθ
=(1)m4π2πg[(1+cosψ)2]{cos[(m+1)ψ]+cos[mψ]}dψ.
(B.5)

That is, bm can be found by averaging adjacent Fourier coefficients from the Fourier series of the periodic and even function 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

9. W. Press, S. Teukolsky, W. Vetterling, and B. Flannery, Numerical Recipes: The Art of Scientific Computing (Cambridge University Press, 1992) Section 2.9.

].

As a minor variation, it is also possible to discretise the integral in Eq. (B.4) by sampling with a spacing of π/(2N) in one of two obvious ways:

bm(1)m2Nj=0N1{11+δj0g[cos2(π2Nj)]cos(π2Nj)}cos[π2N(2m+1)j],
(B.6)
bm(1)m2Nj=0N1{g[cos2(π4N(2j+1))]cos(π4N(2j+1))}cos[π4N(2j+1)(2m+1)].
(B.7)

References and links

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]

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 5494, 113–121 (2004), doi:10.1117/12.551420. [CrossRef]

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]

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, Handbook of Mathematical Functions (Dover, 1978), Chap. 22.

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.

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]

9.

W. Press, S. Teukolsky, W. Vetterling, and B. Flannery, Numerical Recipes: The Art of Scientific Computing (Cambridge University Press, 1992) Section 2.9.

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:  Author  |  Year  |  Journal  |  Reset  

References

  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]
  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 5494, 113-121 (2004), doi:10.1117/12.551420. [CrossRef]
  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]
  5. E. W. Weisstein, "Jacobi Polynomial" from MathWorld—A Wolfram Web Resource. http://mathworld.wolfram.com/JacobiPolynomial.html, see esp. Eqs. (10)-(14).
  6. M. Abramowitz, and I. Stegun, Handbook of Mathematical Functions (Dover, 1978), Chap. 22.
  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.
  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]
  9. W. Press, S. Teukolsky, W. Vetterling, and B. Flannery, Numerical Recipes: The Art of Scientific Computing (Cambridge University Press, 1992) Section 2.9.
  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.

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.

CrossCheck Deposited