OSA's Digital Library

Optics Express

Optics Express

  • Editor: C. Martijn de Sterke
  • Vol. 20, Iss. 3 — Jan. 30, 2012
  • pp: 2483–2499
« Show journal navigation

Characterizing the shape of freeform optics

G. W. Forbes  »View Author Affiliations


Optics Express, Vol. 20, Issue 3, pp. 2483-2499 (2012)
http://dx.doi.org/10.1364/OE.20.002483


View Full Text Article

Acrobat PDF (3456 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

A recently introduced method for characterizing the shape of rotationally symmetric aspheres is generalized here for application to a wide class of freeform optics. New sets of orthogonal polynomials are introduced along with robust and efficient algorithms for computing the surface shape as well as its derivatives of any order. By construction, the associated characterization offers a rough interpretation of shape at a glance and it facilitates a range of estimates of manufacturability.

© 2012 OSA

1. Introduction

Progressive ophthalmic lenses were conceived over 75 years ago [1

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

] and soon became more complex [2

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

,3

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

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

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]

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

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

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]

]. Other interesting design methods are also evolving to enable the discovery of more effective systems involving freeform optics [7

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]

9

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

].

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

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

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]

]. 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=ccon[(x+s)2+y2]1+1(1+κ)ccon2[(x+s)2+y2]=ccon(ρ2+2sρcosθ+s2)1+1(1+κ)ccon2(ρ2+2sρcosθ+s2),
(2.1)
where ccon 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 zaxis 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.
.

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
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.
. Orthogonalization is necessarily performed over a specific domain and, for that purpose, I adopt this circular region that encloses the part.

Notice that the sum on the first line of Eq. (2.2) has precisely the same form as that used in [14

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

] except that the coefficients and polynomials now carry a superscripted zero. This zero distinguishes them from the non-rotationally symmetric components, namely anm, bnm, and Qnm(u2) 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 XjYk where X=ucosθ and Y=usinθ. [It is useful in this to equate real and imaginary parts after expanding the identity um(cosmθ+isinmθ)=(X+iY)m.] The maximum order of the elements that compose the contribution associated with anm or bnm where m>0 is just m+2n, 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+2nT 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 2n+4T.

Without loss of generality, it is possible to assume in Eq. (2.2) that
a01=b01=0.
(2.3)
When they're non-zero, the contribution from these terms to the normal departure is proportional to u(a01cosθ+b01sinθ)(a01x+b01y). 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
n=0Nan1Qn1(1)=n=0Nbn1Qn1(1)=0,andn=0Nan1Qn1(0)=n=0Nbn1Qn1(0)=0. (2.4a,b)
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. Qnm(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

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

] 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,θ):=u2(1u2)n=0Nan0Qn0(u2)+m=1Mumn=0N[anmcos(mθ)+bnmsin(mθ)]Qnm(u2),
(3.1)
the natural extension of the basis used in [12

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

] is to construct Qnm(v) so that the mean square gradient of this normal departure from the best-fit sphere satisfies
|δ(u,θ)|2=(δu)2+1u2(δθ)2=m,n[(anm)2+(bnm)2].
(3.2)
Here, bn0=0 and angle brackets denote the average over the aperture defined by
g(u,θ):=01ππg(u,θ)w(u)udθdu/01ππw(u)udθdu.
(3.3)
The weight function used in [12

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

] is w(u)=[u2(1u2)]1/2 so Eq. (3.3) becomes
g(u,θ):=1π201ππg(u,θ)[1u2]1/2dθdu.
(3.4)
Equations (3.1), (3.2), and (3.4) completely determine Qnm(v).

The trigonometric functions in Eq. (3.1) mean that the different azimuthal orders (i.e. mθ and mθ for mm) 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, Qn0(v) are precisely the (“Qbfs”) polynomials introduced in [13

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

], but those for m>0 require further treatment. After performing the angle integral, the condition for orthonormality reduces to the requirement that
1π01{u[umQnm(u2)],mum1Qnm(u2)}{u[umQnm(u2)],mum1Qnm(u2)}du1u2
(3.5)
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
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,...
. 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
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.
. The complementary basis members, i.e. umsin(mθ)Qnm(u2), have similar plots, of course, but each is rotated by π/(2m).

The plots in Fig. 4 for m=0 follow upon revolving the familiar plots in Fig. 1 of [14

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

], 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 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
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.
, 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×(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].

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:
anmcos(mθ)+bnmsin(mθ)=αnmcos(mθϕnm).
(3.6)
With this, the mean square gradient of Eq. (3.2) satisfies
|δ(u,θ)|2=m,n[(anm)2+(bnm)2]=m,n(αnm)2.
(3.7)
In this sense, αnm 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

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

], 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/cparab=20mm 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=20mm. When the radius of the cylinder is ρmax=10mm 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
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, bnm 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.
and 7
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.
. Notice that symmetry in y means that bnm0. These results use some additional terms beyond those given in Fig. 3:

Q06(x)=897,Q16(x)=891397(7772x), (4.1a,b)
Q07(x)=327231,Q08(x)=8858. (4.1c,d)

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
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, a01 now vanishes. The inverse of the best-fit curvature for this case is 37.432729 mm.
. 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.

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

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

]. In particular, the sag can be expressed as
z=ccon(ρ2+2sρcosθ+s2)1+1(1+κ)ccon2(ρ2+2sρcosθ+s2)+δ(ρ/ρmax,θ)σ(ρ,θ),
(5.1)
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
σ(ρ,θ):=1(1+κ)ccon2(ρ2+2sρcosθ+s2)1κccon2(ρ2+2sρcosθ+s2).
(5.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

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

], but now imagine that only the upper half of that figure holds the off-axis segment of interest. As described in [12

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

], a map of the fringe density in the resulting interferogram is approximately proportional to |δ(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×N pixel grid, it follows from the discussion in [12

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

] that the rms fringe density in Nyquist units, say γ¯, satisfies
γ¯16Nλm,n[(anm)2+(bnm)2].
(5.3)
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

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

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

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

Pnm(x):=(1)n(2n)!!2(2n1)!!Pn(32,m32)(2x1),
(A.1)

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

Pn+1m(x)=[Anm+Bnmx]Pnm(x)CnmPn1m(x),
(A.2)

where the constants are given explicitly by

Anm:=(2n1)(m+2n2)[4n(m+n2)+(m3)(2m1)]/Dnm,
(A.3a)
Bnm:=2(2n1)(m+2n3)(m+2n2)(m+2n1)/Dnm,
(A.3b)
Cnm:=n(2n3)(m+2n1)(2m+2n3)/Dnm,
(A.3c)

with

Dnm:=(4n21)(m+n2)(m+2n3).
(A.3d)

Just as the auxiliary basis did for the polynomials used in [14

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

], Pnm(x) and Eq. (A.2) underpin a number of efficient processes for working with Qnm(x). The scaling factor in Eq. (A.1) is chosen so that the peak absolute value of umPnm(u2) over 0u1 never exceeds about unity for all m and n. As in [14

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

], this makes it straightforward to determine when terms are negligible and can be dropped from a linear combination.

For any m1, the recurrence relation in Eq. (A.2) would normally be initialized by using P0m(x)=1/2 and P1m(x)=m12(m1)x. [In fact, just the first of these is sufficient for m>3 since C0m=0 and Eq. (A.2) itself generates P1m(x) by taking n=0.] Notice that for m=1, however, P1m(x) would then not be a linear function but a constant. Further, Dnm vanishes when m=n=1 thereby causing Eqs. (A.2) and (A.3) to fail in the generation of P21(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

Pnm(x):={1x/2,m=1andn=1,(1)n(2n)!!2(2n1)!!Pn(32,m32)(2x1),otherwise.
(A.4)

In this case, the first two polynomials are

P0m(x)=12andP1m(x)={1x/2,m=1,m12(m1)x,m>1. (A.5a,b)

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

P21(x)=[3x(128x)]/6andP31(x)={5x[60x(12064x)]}/10. (A.6a,b)

If the integral in Eq. (3.5) is taken to define the inner product written as Qnm,Qnm then it turns out that, for any fixed m1, Pnm,Pnm is a symmetric tri-diagonal matrix. The proof of this follows upon applying five standard identities, namely

x[Pn(a,b)(2x1)]=(n+a+b+1)Pn1(a+1,b+1)(2x1),
(A.8)
Pn(32,m32)(z)=1m+2n2[(m+n2)Pn(12,m32)(z)(m+n32)Pn1(12,m32)(z)],
(A.9)
01xm32[Pn(12,m32)(2x1)]2dx1x=π(2n1)!!(2m+2n3)!!2m+2n1(m+2n1)n!(m+n2)!,
(A.10)
Pn(12,m32)(z)=1m+2n1[(m+n1)Pn(12,m12)(z)+(n12)Pn1(12,m12)(z)],
(A.11)
xPn(12,m12)(2x1)={(n+1)(m+n)(m+2n+1)(m+2n)Pn+1(12,m12)(2x1)+[4mn+m(2m1)+4n21]2(m+2n+1)(m+2n1)Pn(12,m12)(2x1)+(2n1)(2m+2n1)4(m+2n)(m+2n1)Pn1(12,m12)(2x1)}.
(A.12)

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 Pnm,Pnm is zero whenever |nn|>1. The diagonal elements, say Fnm:=Pnm,Pnm, turn out to be given by

Fnm={m2(2m3)!!2m+1(m1)!,n=0,4(n1)2n2+18(2n1)2+1132δn1,n>0andm=1,2nχ(35m+4nχ)+m2(3m+4nχ)(m+2n3)(m+2n2)(m+2n1)(2n1)γnm,n>0andm>1,
(A.13)

where χ=m+n2 as a shorthand, δmn is the Kronecker delta, and

γnm:=n!(2m+2n3)!!2m+1(m+n3)!(2n1)!!.
(A.14)

Similarly, the off-diagonal elements, say Gnm:=Pnm,Pn+1m, are found to be

Gnm={(2m1)!!2m+1(m1)!,n=0,(2n21)(n21)8(4n21)124δn1,n>0andm=1,[2n(m+n1)m](n+1)(2m+2n1)(m+2n2)(m+2n1)(m+2n)(2n+1)γnm,n>0andm>1.
(A.15)

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

γ1m+1=2m+12(m1)γ1mandγn+1m=(n+1)(2m+2n1)(m+n2)(2n+1)γnm.
(A.16)

Of course, a similar process can be used for evaluating F0m and G0m.

Pnm(x)={f0mQ0m(x),n=0,fnmQnm(x)+gn1mQn1m(x),n>0,
(A.17)

where the constants fnm and gnm are found by Cholesky decomposition of Pnm,Pnm. That is, for each azimuthal order —i.e. for any fixed m— these constants are determined by starting with f0m=F0m and working progressively up from n=1 to any desired terminal value by using

gn1m=Gn1m/fn1m,andfnm=Fnmgn1mgn1m. (A.18a,b)

The significance of Eqs. (A.2) and (A.17) is discussed further in Appendix B. For efficiency, Anm, Bnm, Cnm, fnm and gnm 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

[F01F11F21F31F02F12F22F32F03F13F23F33F04F14F24F34]=[141532177229401278353667402732351635162438054511128236122872560],[G01G11G21G01G12G22G01G13G23G01G14G24]=[14124740385487161532732117160356421643332],
(A.19)
[f01f02f03f04f11f12f13f14f21f22f23f24f31f32f33f34]=[121234325214721419214185634273216115141214538141280337012278518315339723014684129094141132560611612890571114],
(A.20)
[g01g02g03g04g11g12g13g14g21g22g23g24]=[123425467532131456387433701276121107230741929011743712803033161832785].
(A.21)

The results in this appendix give ready access to Qnm(x) at any order. In particular, it follows from Eq. (A.17) that we can start with Q0m(x)=P0m(x)/f0m=1/(2f0m) and then work up in n by using

Qnm(x)=[Pnm(x)gn1mQn1m(x)]/fnm.
(A.22)

Keep in mind that Pnm(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 Qnm(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

Pnm(0)={1,m=1andn=1,(2m+2n3)!!2(2m3)!!(2n1)!!,otherwise.
(A.23)

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

f(ρ,θ)=cρ21+1c2ρ2+11c2ρ2{u2(1u2)n=0Nan0Qn0(u2)+m=1Mum[cos(mθ)n=0NanmQnm(u2)+sin(mθ)n=0NbnmQnm(u2)]}.
(B.1)

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

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

]. For any fixed m1, this appendix addresses the additional core task of working with the sums of the form of those on the second line, say

S(x):=n=0NcnQnm(x),
(B.2)

where cn may be either anm or bnm. (To reduce clutter, neither S nor cn takes a superscripted m.) By following the same procedure discussed in Section 3 of [14

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

], it follows from Eq. (A.17) above that S(x) can also be expressed in terms of the auxiliary polynomials as

S(x)n=0NdnPnm(x),
(B.3)

where dn can be found simply from cn by starting with dN=cN/fNm and working down in n from N1 to 0 with

dn=(cngnmdn+1)/fnm.
(B.4)

The inverse of this process is also sometimes useful: for n from 0 to N1 use

cn=fnmdn+gnmdn+1
(B.5)

and finish with cN=fNmdN. That is, for any m>0, superpositions of Qnm(x) and Pnm(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

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

] are directly applicable. In particular, S(x) can be evaluated by using Clenshaw’s method, which starts with αN=dN and αN1=dN1+(AN1m+BN1mx)αN and progressively works down from n=N2 with

αn=dn+(Anm+Bnmx)αn+1Cn+1mαn+2.
(B.6)

(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)=12α0.

When m=1, on the other hand, it helps to choose values for An1, Bn1, and Cn1, 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

A01=2,B01=1,A11=4/3,B11=8/3,C11=11/3,C21=0,
(B.7)
A0m=2m1,B0m=2(1m),whenm>1,
(B.8)

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

S(x)={12α025α3,m=1andN>2,12α0otherwise.
(B.9)

Notice that αn of Eq. (B.6) is a polynomial in x of order Nn.

One of the most striking aspects of Eqs. (B.6) and (B.9) is that, as described in Appendix B of [16

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

], 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)={12α0(j)25α3(j),m=1andN>2,12α0(j)otherwise,
(B.10)
αn(j)=jBnmαn+1(j1)+(Anm+Bnmx)αn+1(j)Cn+1mαn+2(j),whenj>0.
(B.11)

When initialized with αN+2j(j)=αN+1j(j)=0, Eq. (B.11) can proceed downward from n=Nj to give the desired results for Eq. (B.10). Notice that this process for evaluating the j’th derivative assumes that the (j1)’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 cn to dn 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

[A01A02A03A04A05A06A11A12A13A14A15A16A21A22A23A24A25A26A31A32A33A34A35A36A41A42A43A44A45A46]=[235791143232762785241062595261521175020375108355528663525362451355613049161811226325152431436316166],
(B.12a)
[B01B02B03B04B05B06B11B12B13B14B15B16B21B22B23B24B25B26B31B32B33B34B35B36B41B42B43B44B45B46]=[12468108344409528524544215112252453074414435307220491122744110278821133],
(B.12b)
[C01C02C03C04C05C06C11C12C13C14C15C16C21C22C23C24C25C26C31C32C33C34C35C36C41C42C43C44C45C46]=[******1133533527987775059715215088225133527282125273589112253956334980814549556314301701404911051386].
(B.12c)

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]

OCIS Codes
(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: November 21, 2011
Revised Manuscript: January 4, 2012
Manuscript Accepted: January 6, 2012
Published: January 19, 2012

Citation
G. W. Forbes, "Characterizing the shape of freeform optics," Opt. Express 20, 2483-2499 (2012)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-20-3-2483


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  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. SPIE8167, 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. SPIE7061, 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. SPIE7652, 76520C, 76520C-7 (2010). [CrossRef]
  8. J. R. Rogers, “A comparison of anamorphic, keystone, and Zernike surface types for aberration correction,” Proc. SPIE7652, 76520B, 76520B-8 (2010). [CrossRef]
  9. A. Yabe, “Method to allocate freeform surfaces in axially asymmetric optical systems,” Proc. SPIE8167, 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. SPIE8169, 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. Express19(10), 9923–9941 (2011). [CrossRef] [PubMed]
  13. G. W. Forbes, “Shape specification for axially symmetric optical surfaces,” Opt. Express15(8), 5218–5226 (2007). [CrossRef] [PubMed]
  14. G. W. Forbes, “Robust, efficient computational methods for axially symmetric optical aspheres,” Opt. Express18(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. Express15(26), 18014–18024 (2007). [CrossRef] [PubMed]
  16. G. W. Forbes, “Robust and fast computation for the polynomials of optics,” Opt. Express18(13), 13851–13862 (2010). [CrossRef] [PubMed]

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