Open Access
26 January 2016 Aspherical surface profile fitting based on the relationship between polynomial and inner products
Author Affiliations +
Abstract
High-precision aspherical polynomial fitting is essential to image quality evaluation in optical design and optimization. However, conventional fitting methods cannot reach optimal fitting precision and may somehow induce numerical ill-conditioning, such as excessively high coefficients. For this reason, a projection from polynomial equations to vector space was here proposed such that polynomial solutions could be obtained based on matrix and vector operation, so avoiding the problem of excessive coefficients. The Newton–Raphson iteration method was used to search for optimal fitting of the spherical surface. The profile fitting test showed that the proposed approach was able to obtain results with high precision and small value, which solved the numerical ill-conditioning phenomenon effectively.

1.

Introduction

With the rapid development of optical design and manufacturing technology, aspherical profiles have seen extensive use because of their lower rate of aberrations relative to spherical profiles.14 Improving the precision of fitting projects would facilitate the evaluation of image quality in optical design and optimization significantly. This is commonly done by increasing the number of samples or using high-order polynomials. Currently, there are a number of articles that discuss improving the precision of polynomial fitting. Because increasing the number of samples beyond a certain point will not markedly improve fitting precision,5 high-order polynomials have seen widespread use. The least square method will usually produce ill-conditioned Gram matrixes and so introduce large errors in solution. To solve the problem, the Gram–Schmidt process was used to generate mutually orthogonal polynomials in the unit circle.6,7 The Gauss–Newton algorithm was used to solve nonlinear least squares problems through iterations to reduce fitting errors. This method was reported to be able to further improve fitting precision.8,9

However, when fitting rotational symmetric surfaces using the nonorthogonal polynomials and least square method, the peak polynomial coefficient increases with increasing polynomial order such that massive coefficients may be derived in some cases.10 These coefficients are usually randomly oriented and several orders of magnitude higher than rise of arch, which might result in a much smaller number from the subtraction of two large numbers when calculating the rise of arch. In this way, the limited accuracy of numerical storage in computers leads to losses of effective numbers during computation. Specifically, it compromises the accuracy of the calculation of the surface slope, thereby rendering ray tracing inaccurate. In optical optimization processes, the loss of effective numbers causes numeric instability of the optimization algorithm, which is unacceptable in practical applications.11 Although the methods listed above could improve fitting precision, it is not possible to prevent the loss of effective numbers because of large coefficients. In this way, a projection from polynomial equations to vector space was proposed in this study, thereby introducing a new fitting method, achieving the optimal quadric fitting surface with better fitting errors as well as minimizing the magnitude of the aspheric coefficient.

In order to improve fitting precision effectively and solve the problem of excessive coefficients, it is necessary to analyze the sources of fitting errors and their relationship with the coefficients of polynomials. The coefficient of each order of polynomials is not directly correlated to fitting errors, so direct analysis is not efficient. Based on the characteristics of homogeneity and additivity of polynomials and vectors and on the similar definition of inner products and the simple linear relationship between components and modulus length of each vector and error in least square calculation, the problem can be solved by projecting polynomials into vectors. It is inevitable to calculate the inner products of polynomials while conducting profile fitting with either the linear or the nonlinear least square method. Because the definition of fitting error is directly related to inner products, the projection from polynomials to vector space can be built based on its relationship with inner products such that the polynomials are analyzed using vector analysis methods, thereby transforming the polynomial problem into a vector problem.

2.

Methods

2.1.

Projection of Polynomials into Vector Space

The polynomial equations consisting of n polynomials were here projected into n-dimension vector space, with each polynomial projected into the n-dimension vector. The polynomials are written as ϕ1,ϕ2,,ϕn, and the corresponding vectors as v1,v2,,vn. In order to facilitate future calculations, the following projection principles were here adopted:

Eq. (1)

(1)(φi,φi)=(vi,vi),i=1,2,,n(2)(φi,φj)=(vi,vj),i,j=1,2,,n,ij(3)vi=(mi1,mi2,,mii,0,,0)T.
The vectors are integrated as matrix M

Eq. (2)

M=[v1,v2,,vn].

One algorithm is to obtain matrix M that satisfies the projection principles through “Cholesky” decomposition of the “Gram” matrix in Eq. (3). From the Cholesky decomposition, a lower triangular matrix and an upper triangular matrix, which are mutual symmetric, could be obtained, where M is the upper triangular matrix. The polynomials coefficient matrix derived from fitting is recorded as A, and the projected vector is v. The coefficient matrix and the corresponding vector have the relationships as

Eq. (3)

v=MA

Eq. (4)

A=M1v,
where A=(A1,A2,,An)T.

2.2.

Principles of Aspheric Polynomial Surface Fitting

In aspheric polynomial surface fitting, the profile expression conforming to ISO 10110 standard is adopted12

Eq. (5)

z=cr21+1(1+k)c2r2+j=0Ajr2j+4.
Here, r=x2+y2 is the perpendicular distance from surface points to the optic axis; z is the rise of the arch, i.e., the perpendicular distance between surface points to the xy plane; c is the curvature of surface apex; k is the conical constant; and Aj is the polynomial coefficient.

Using optimal spherical surface methods, traditional fitting methods, based on surface features, first employed apex curvature and edge rise of arch to calculate quadric surface parameters c and k.13,14 d(r): d(r)=z(r)cr21+1(1+k)c2r2 is the difference between fitting surface z(r) and optimal fitting quadric surface.

The least square method is also generally used to fit residuals. The polynomials are represented as ϕ0(r), ϕ1(r),,ϕn1(r) in rectangular coordinate systems; the fitted residual function is f(r), and W(r) is a weight function. The least square calculation is used to derive A0,A2,,An1, which makes (dAjr2j+2)2 minimal. The solution equations are as follows:9,10

Eq. (6)

GA=b.
The coefficients of every polynomial A1,A2,,An, are obtained by solving Eq. (6), where G is the Gram matrix and its elements are given by Eq. (7), and b is the calculation matrix and defined by b=[d·r4,,d·r2n+2]T. And A is the vector of coefficients defined by A=[A1,,An]T.

Eq. (7)

Gij=r2i+4r2j+4.

2.3.

Optimal Quadric Fitting Surface

As stated in Sec. 2.2, during the fitting process, the optimal quadric fitting surface is determined first, and then the least square method is used to fit the residuals. Using polynomials with more terms involves more dimensions of the corresponding vector space. The data to be fitted can be said to correspond to an infinite dimensional vector. The optimal quadric fitting surface is also represented by an infinite dimensional vector, as shown in Eq. (8), where vf is the corresponding vector of the surface to be fitted, and vs is the vector of the optimal quadric fitting surface.

Eq. (8)

vf=limi[α0,α1,,αn1,αn,,αi]Tvs=limi[β0(c,k),β1(c,k),,βn1(c,k),βn(c,k),,βi(c,k)]T,
where αi and βi are the values of vector components, and βi represents the component values of the vector of the optimal quadric fitting surface, which could be represented as the function of c and k. The vector vd corresponding to residuals is shown

Eq. (9)

vd=limi(α0β0,α1β1,,αiβi)T.

During the fitting process, the use of one additional polynomial produces an additional vector component, and the fitting error ve is proportional to the component modulus value in the undescribed dimension

Eq. (10)

ve=limi(0,,αnβn,αn+1βn+1,,αiβi)T.

When the components of vector ve that exceeds n dimensions have minimal modulus value, the modulus length of ve is minimal, and the fitting precision is maximal. In this way, it can be deduced that there is an optimal quadric fitting surface that can yield the minimal modulus value of undescribed components

The determination of optimal quadric fitting surface is a nonlinear least square problem, which is commonly solved using the Gauss–Newton algorithm. The principle is to perform the first-order Taylor expansion on the base function to linearize it. Then the linear least square method is used in fitting, and the final results can be obtained upon iteration convergence. However, this method has several drawbacks. For example, it uses local quadratic convergence for the zero residual fitting problem, the convergence speed is low for small residual fitting problems, and there is no convergence for large residual fitting problems. Here, the residual refers to the discarded high-order terms in Taylor expansion. For this reason, the Newton–Raphson method was used in this study to perform the second-order Taylor expansion based on the definition formula of fitting residual. Because one more order of Taylor expansion was performed than in the Gauss–Newton method, more information from the function was preserved, so the method showed good convergence. The procedure is as follows:

The square sum s of fitting residual d is defined as follows: d=zcr21+1(1+k)c2r2i=0Ajr2j+4, s=d2. The target of the algorithm is to find c and k that conform to sc=sk=0. First, the initial values of c and k are obtained based on conventional methods. Then, Eq. (11) was used in iteration so as to obtain new values of c and k.

Eq. (11)

[cn+1kn+1]=[cnkn][2sc22sck2sck2sk2]1[scsk].

Because in the fitting process, Aj is determined once c and k are determined, Aj is correlated to c and k.

Eq. (12)

sc=2d(dcaicr2i+4),[a0can1c]=[(φ0,φ0)(φ0,φn1)(φn1,φ0)(φn1,φn1)]1[(dc,r4)(dc,r2n+2)].

The partial derivative of k is obtained with similar method. The second-order partial derivative matrix was obtained based on the calculation of the first-order partial derivative. New values of c and k are continuously obtained according to the iteration formula until the iteration termination standard is achieved.

2.4.

Methods to Decrease Polynomial Coefficients

Once the optimal quadric fitting surface is determined, the least square method is used to find the polynomial coefficients. If the polynomials have a large number of terms or the surface space frequency is high, there might be too many coefficients in the fitting results, and these can be several orders of magnitude higher than rise of arch. This is unfavorable in practical applications, so it is necessary to decrease the polynomials coefficients without affecting the fitting precision. The solving principle of the least square method involves determining the minimal root-mean-square (RMS) error, so the obtained small fitting coefficients can cause increased RMS (the least square method assures minimal RMS, but the coefficients can be large, and none of the operations on the coefficients assure minimal RMS). As shown, the high fitting precision and small coefficients cannot be obtained simultaneously. In this way, small coefficients could be obtained when a low fitting precision is acceptable. Proper coefficients could be obtained by limiting RMS and the range of polynomial coefficients. However, solving the multivariable quadratic inequality is difficult and time consuming. The size of the coefficients could be reduced by controlling the variation in polynomial coefficients and error-associated variables.

As stated above, additive error can be introduced while reducing polynomial coefficients. The ΔA change in the polynomial coefficients matrix introduces error E:

Eq. (13)

E=ΔAifi2=ΔAivi2=MΔA2.
Here, ΔA=(ΔA1,ΔA2ΔAn)T, and ΔAi is the coefficient variation of the i polynomial. It is known from Eq. (4) that ΔA=M1Δv. In general, the absolute values of M1 row elements increases from the left to the right so the last component of Δv makes the largest contribution to the coefficient. In other words, given the same error, it can, to the greatest extent, change the polynomial coefficients if the error is allocated to the last component of the vector.

For this reason, the unit error vector is written as u=(001)T, and the variation of polynomial coefficients is B=kM1u when k times the unit error is introduced. Thus, the coefficient of new polynomials is A+B. The introduced error can be determined by taking minimal square sum of polynomials coefficients as the standard for coefficient reduction

Eq. (14)

k=(A,B)(B,B)=(A,M1u)(M-1u,M1u).

Finally, the coefficient matrix of new polynomials can be generated

Eq. (15)

A=A+kM1u=A(A,M1u)(M1u,M1u)M1u.

The fitting error of new polynomials can be calculated based on additive error and fitting error of original polynomials other than sample points’ coordinates. From the perspective of vector dimensions, the fitting error of original polynomials is due to finite terms. For this reason, the dimension of corresponding vector is finite, resulting in undescribed high-dimension features. It is known that different dimensions of the corresponding vector of the additive error and original fitting error are orthogonal, which allows the errors to be integrated based on mean square synthesis method.

3.

Examples

In this section, the examples are illustrated. The aplanatic Descartes oval profile was used to validate the feasibility, discussing the departure part, and the optimal quadric fitting surface part. Another two examples, an extreme and a mild aspheric profile, will perform comprehensive validation of the method. In the calculation procedure, the Eq. (15) was applied after the optimal quadric fitting surface was determined by using the proposed method and the polynomial coefficients were solved with the least square method.

3.1.

Aplanatic Descartes Oval Profile

The aplanatic Descartes oval profile which used by Forbes was used here to validate the feasibility of the method.10 Using the conventional method, standard aspheric polynomial fitting was performed. Results are shown in Table 1. The curves of Descartes oval and their fitting errors are shown in Fig. 1. As stated in Forbes’s article, standard aspheric polynomials are nonorthogonal, and the use of high-order terms in the fitting process can affect the coefficients of low-order terms, resulting in higher coefficients from more terms. In this way, the number of significant figures can decrease during fitting with more high-order polynomials, leading to nonsignificant increases in precision. In this paper, further analysis and deduction were conducted, and we demonstrate that the results were not optimal either in terms of precision or in terms of coefficient magnitude. The effectiveness of our proposed method was proved using two conditions: (1) excessive coefficients and (2) suboptimal fitting precision from optimal quadric fitting surface.

Table 1

Fitting result using aspheric polynomials for Descartes oval.

No. of termsTraditional methodImprovement for departureImprovement for optimal spherical fitting surface
FiveSevenEightEightEight
c (mm1)0.50.50.50.50.5
k0.8437160.8437160.8437160.8437160.759214
A1 (mm)2.57757×1012.69713×1012.73894×1012.71357×1014.70730×101
A2 (mm)5.3405×1021.45314×1012.19449×1011.74068×1017.2462×102
A3 (mm)4.33431×1012.72898×1012.28076×1018.1338×1028.9275×102
A4 (mm)5.73476×1018.36713×1018.81902×1011.89150×1014.788×103
A5 (mm)3.44232×1011.45812×1001.834551×1002.36148×1012.7192×102
A6 (mm)1.19517×1012.374638×1001.09050×1017.1639×102
A7 (mm)4.25257×1011.622536×1003.10880×1015.6685×102
A8 (mm)4.83067×1011.70768×1012.2460×102
RMS(mm)4.36995×1053.5701×1069.8508×1072.6815×1061.2038×108

Fig. 1

Descartes oval (pink line “a”: Descartes oval, green line “b”: the error of improvement method for departure, blue line “c”: the error of traditional method, red line “d”: the error of improvement method for optimal spherical fitting surface).

OE_55_1_015105_f001.png

The data in the table show that, when there are eight polynomial terms, the coefficients are large, two orders of magnitude greater than the maximal value of the point of departure. These very large coefficients affect the system optimization process and so need to be decreased. According to the proposed method, taking the eight-term aspheric polynomials, e.g., the polynomials were projected into eight-dimensional vector space, and the eight vectors v1,v2v8 are here represented as matrix M

M=[4.071e+03.378e+02.898e+02.546e+02.278e+02.066e+01.895e+01.754e+006.239e19.275e11.078e+01.150e+01.179e+01.185e+01.176e+0001.245e12.590e13.747e14.675e15.398e15.953e10002.733e27.217e21.242e11.769e12.269e100006.271e31.997e23.969e26.345e2000001.474e35.478e31.232e20000003.513e41.489e300000008.435e5].

The associated variables between polynomial coefficients variation and introduced error can be defined and controlled for optimization, facilitating reduction of coefficient size. The results showed that the coefficients decreased into the same level as in fitting of five-term polynomials, while the fitting precision was kept at a level similar to that of the fitting of seven-term polynomials.

Similarly, based on the deduction shown in Sec. 3, the Newton–Raphson iteration method was used to search for optimal fitting spherical surface. Results are shown in Table 1. The optimal fitting spherical surface is distinct from that produced using conventional methods. After determining the optimal quadric fitting surface using the Raphson method, the polynomial coefficients were solved, and the fitting error was dramatically decreased.

3.2.

Comprehensive Validation Using an Aspheric Surface

The examples given above are experimental validation of the new method in terms of the departure part and the optimal quadric fitting surface part. The following example will perform comprehensive validation for the above two parts. The surface as shown in Fig. 2 is the measured contour from a fabricated lens used in a cell phone. It was fitted with the proposed method, and the results are as shown in Table 2. The results demonstrated that, due to the variation in the optimal quadric fitting surface, the vector corresponding to the point of departure shows high-dimensional components with lower values, resulting in more described surface information in the polynomials, so retaining high fitting precision. Associated variables between polynomials coefficients variation and introduced error were obtained, and this could significantly change coefficient range while keeping the error similar. The results also show that the proposed method is effective for a measured profile as well as the nominal starting values of a surface.

Fig. 2

An aspheric surface profile used in cellphone lens (pink line “a”: the aspheric surface, blue line “b”: the error of traditional method, green line “c”: the error of the proposed method).

OE_55_1_015105_f002.png

Table 2

Fitting result using aspheric polynomials for an aspheric surface.

Number of termsTraditional methodImproved method
FiveEightEight
c (mm1)6.644046.644044.83301
k2.42266×1012.42266×1011.50273×100
A1 (mm)1.22806×1012.89211×1011.66517×101
A2 (mm)4.92324×1012.37090×1027.30280×100
A3 (mm)7.79851×1019.09388×1022.00049×101
A4 (mm)5.59258×1011.97371×1032.01305×101
A5 (mm)1.50155×1012.54711×1033.30800×102
A6 (mm)1.93555×1031.56877×101
A7 (mm)7.99407×1021.16211×101
A8 (mm)1.38363×1022.74648×100
RMS (mm)2.39863×1027.51581×1036.21144×104

3.3.

Validation Using a Mild Aspheric Surface

The examples in Secs. 3.1 and 3.2 are extreme aspheric surfaces. In this section, the proposed method is also applied to the aspheric surface fitting of a mild profile, the contour from a imaging lens. The fitting results were shown in Fig. 3 and Table 3, which approved the efficiency of the method.

Fig. 3

The aspheric surface of a mild profile (red line “a”: the aspheric surface, blue line “b”: the error of traditional method, green line “c”: the error of the proposed method).

OE_55_1_015105_f003.png

Table 3

Fitting result using aspheric polynomials for an aspheric surface of mild profile.

Number of termsTraditional methodImproved method
FiveEightEight
c (mm1)1.000931.000930.737635
k2.18644×1002.18644×1006.60064×101
A1 (mm)1.37466×1003.11691×1001.09348×100
A2 (mm)4.45997×1002.19149×1013.05223×100
A3 (mm)5.95359×1007.45048×1014.91347×100
A4 (mm)3.66472×1001.43932×1023.33001×100
A5 (mm)8.48571×1011.65021×1021.09858×100
A6 (mm)1.11146×1023.27368×100
A7 (mm)4.06120×1011.94461×100
A8 (mm)6.21157×1003.93294×101
RMS (mm)3.22914×1031.20238×1031.874801×104

4.

Conclusions

While performing surface fitting with the least square method, polynomial coefficients often become far larger than maximal rise of arch of the fitted data. This causes difficulties in optical system design and optimization. In this paper, a projection from polynomials to vector was proposed such that the polynomial problem was transformed into a vector problem. According to the proposed projection principle, the mutual projection between polynomials and vector could determine the relationship between the vector and the polynomial coefficients, and further determine the relationship between introduced errors and polynomial coefficient variation. The polynomial coefficients were reduced based on the standard of minimal square sum. For the first time, the Newton–Raphson method was used to search for optimal fitting surface, which not only improved the fitting precision but also had better convergence than the Gauss–Newton method. Through fitting tests of various surface profiles, the proposed method was proven to be effective in terms of fitting error and the range of the coefficients. Thus, it is of significance in practical application. The proposed algorithm has been programmed with MATLAB and available for a Dynamic Data Exchange connection to optical design software.

Acknowledgments

This work was supported in part by the grant from the National Natural Science Foundation of China (Nos. 61275003, 51327005, and 91420203).

References

1. 

A. Yoshikazu et al., “On-machine measurement of aspherical surface profile,” Nanotechnol. Precis. Eng., 2 (3), 210 –216 (2004). Google Scholar

2. 

D. Huang et al., “Shape measurement of aspheric plastic lens with large angle,” Chin. Opt. Lett., 3 (9), 516 –519 (2005). CJOEE3 1671-7694 Google Scholar

3. 

R. Q. Wu et al., “Analysis of athermalizing performance of thermal infrared with Cassegrain antenna,” Optik, 121 1904 –1907 (2010). http://dx.doi.org/10.1016/j.ijleo.2009.05.012 OTIKAJ 0030-4026 Google Scholar

4. 

S. J. Ludwick et al., “Design of a rotary fast tool servo for ophthalmic lens fabrication,” Prec. Eng., 23 253 –259 (1999). http://dx.doi.org/10.1016/S0141-6359(99)00017-3 Google Scholar

5. 

S. L. Xie and Y. Zhong, “Research of sampling point number in wave front fitting with Zernike polynomials,” Laser Technol., 35 (2), 272 –277 (2011). Google Scholar

6. 

D. Hilbig et al., “Fitting discrete aspherical surface sag data using orthonormal polynomials,” Opt. Exp., 23 (17), 22404 –22413 (2015). http://dx.doi.org/10.1364/OE.23.022404 Google Scholar

7. 

G. W. Forbes, “Characterizing the shape of freeform optics,” Opt. Exp., 20 (3), 2483 –2499 (2012). http://dx.doi.org/10.1364/OE.20.002483 OPEXFF 1094-4087 Google Scholar

8. 

W. Sun, J.W. McBride and M. Hill, “A new approach to characterising aspheric surfaces,” Prec. Eng., 34 171 –179 (2010). http://dx.doi.org/10.1016/j.precisioneng.2009.05.005 Google Scholar

9. 

M. L. N. Gonçalvesa and P. R. Oliveirab, “Convergence of the Gauss–Newton method for a special class of systems of equations under a majorant condition,” Optimization, 64 (3), 577 –594 (2015). http://dx.doi.org/10.1080/02331934.2013.778854 OPTZDQ Google Scholar

10. 

G. W. Forbes, “Shape specification for axially symmetric optical surfaces,” Opt. Exp., 15 (8), 5218 –5226 (2007). http://dx.doi.org/10.1364/OE.15.005218 OPEXFF 1094-4087 Google Scholar

11. 

G. W. Forbes, “Robust and fast computation for the polynomials of optics,” Opt. Exp., 18 (13), 13851 –13862 (2010). http://dx.doi.org/10.1364/OE.18.013851 OPEXFF 1094-4087 Google Scholar

12. 

G. H. Spencer and M. V. R. K. Murty, “General ray-tracing procedure,” J. Opt. Soc. Am., 52 (6), 672 –678 (1962). http://dx.doi.org/10.1364/JOSA.52.000672 JOSAAH 0030-3941 Google Scholar

13. 

H. Aceves-Campos, “Profile identification of aspheric lenses,” Appl. Opt., 37 (34), 8149 –8150 (1998). http://dx.doi.org/10.1364/AO.37.008149 APOPAI 0003-6935 Google Scholar

14. 

Z. Zhang, “Parameter estimation techniques: a tutorial with application to conic fitting,” Image Vision Comput., 15 59 –76 (1997). http://dx.doi.org/10.1016/S0262-8856(96)01112-2 Google Scholar

Biography

Xuemin Cheng is an associate professor at the Graduate School at Shenzhen, Tsinghua University. She received her PhD in optical engineering from Beijing Institute of Technology, China, in 2004. She is author of 56 peer reviewed publications related to optical design, vision, and image quality assessment. Her current research interests include optical design and engineering, including Imaging system design, illumination system design, optical modeling and simulation, etc. She is a member of SPIE.

Yikang Yang received his bachelor’s degree in precise instrumentation from Tianjin University, Tianjin, China, in 2013. He is currently a master's degree candidate at the Tsinghua University.

Qun Hao obtained her PhD in optical instruments from Tsinghua University, Beijing, China in 1998. She is currently a professor in the Department of Optoelectronics of Beijing Institute of Technology (BIT) in Beijing, China, and served as dean of the School of Optoelectronics at BIT.

CC BY: © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.
Xuemin Cheng, Yikang Yang, and Qun Hao "Aspherical surface profile fitting based on the relationship between polynomial and inner products," Optical Engineering 55(1), 015105 (26 January 2016). https://doi.org/10.1117/1.OE.55.1.015105
Published: 26 January 2016
Lens.org Logo
CITATIONS
Cited by 4 scholarly publications.
Advertisement
Advertisement
KEYWORDS
Aspheric lenses

Error analysis

Vector spaces

Spherical lenses

Optical design

Optical engineering

Chromium

RELATED CONTENT

Direct design of a two surface lens including an entrance...
Proceedings of SPIE (September 23 2015)
Optical system design and measurement of freeform surface
Proceedings of SPIE (October 10 2020)
New method for calculating third , fifth , and seventh...
Proceedings of SPIE (December 01 1991)
Nonlinearity and lens design
Proceedings of SPIE (January 01 1991)

Back to Top