- Email: [email protected]
Vibration monitoring by fiber optic fringe projection and Fourier transform analysis
15 June 1997 OPTICS COMMUNICATIONS ELSEVIER
Optics Communications I39 (1997) 17-23
Vibration monitoring by fiber optic fringe projection and Fourier transform analysis Giuseppe Schirripa Spagnolo, Domenica Paoletti, Dario Ambrosini llipurtimento
di Energetica.
Uniwrsifir
drgli
Studi di L’Aquila.
Locali&
Mont&co
di Rob,
67040 Roio Poggio
- L’Aquila,
Italy
Received 14 November 1996; revised 15 January 1997; accepted 16 January 1997
Abstract A fiber optic fringe projection interferometer is proposed to obtain local amplitude vibration of a diffuse surface. A sinusoidal fringe pattern is projected onto the vibrating surface by an optical fiber version of Young’s interferometer. An image of this surface is captured by a CCD camera and fed to a PC through a frame grabber and hence displayed onto a TV monitor. The visibility of the interference pattern is modulated by a function of the local amplitude of vibration. Thus loci of constant vibrational amplitude can be directly observed as regions of low fringe visibility in the image of the surface. A quantitative analysis is performed by a Fast Fourier Transform (FFT) algorithm. Experimental results in agreement with the theory are presented. Keywords:
Vibration analysis; Fiber optic; Digital image analysis: Interferometry; Non-destructive testing
1. Introduction Many coherent optical techniques, based on holographic and speckle interferometry and its variants, have been devised for experimentally studying the vibrating surfaces [l-5]. Recently the introduction of digital image processing methods has led to new possibilities in this field, especially in the application of the speckle interferometry, but, despite any technical improvements, the limits to the quality of speckle images and the complex technologies of the holographic systems (the stringent stability requirements, etc.) are significant impediments to realize systems suitable for routine measurements in industrial environments. In this paper we propose to study vibrating objects by a sinusoidal fringe projection interferometer just developed for contour maps [6-91. The system is based on a periodic fringe pattern projected onto a vibrating surface which is imaged by a zoom lens onto the photosensor of a CCD camera. In most cases, the vibration frequency is much higher than the framing frequency of the CCD camera, so that the video signal, displayed on a monitor. is a time-averaged intensity pat-
tern. The resulting fringe pattern on the vibrating surface carries information regarding the local vibrational amplitude. Loci of constant vibrational amplitude can be directly observed as regions of high or low fringe visibility in the image of the surface. To get quantitative information on the local amplitude of vibration a Fast Fourier Transform (FIT) algorithm can be performed on the time-averaged CCD image. The FFI provides a high-quality display of J,, fringes that can be used to determine the vibration amplitude. The system is easy and cheap, with low mechanical stability requirements. It is well suited to measure vibrations with relatively large amplitude. The proposed method can be considered to some extent a new electronic version of the technique described by Vest and Sweeney [lo].
2. Experimental
setup and theory
The fringe projection system proposed in this paper is based on an optical fiber version of the well-known
0030-4018/97/$17.00 Copyright 0 1997 Elsevier Science B.V. All rights reserved. PII soo30-40 I X(97)00045-X
G. Schirripcr Spngnolo et al. / Optics Communications
Young’s interferometer [I I]. Interference fringes are formed with the coherent light of two single-mode fibers, which replace the pair of pin-holes in the Young’s configuration. In the classical version the power loss associated with the pin-holes is compensated by focusing the light onto their apertures. By doing so, however, Airy diffraction disks become visible and the regularity of the fringe pattern is altered as a consequence. These disadvantages are eliminated by the introduction of optical fibers, which allow Young’s fringes to be projected with a low power laser source.
139 11997) 17-23
The schematic diagram of the system is given in Fig. I. A graded-index (grin) lens is used to couple the light of a 3.5 mW HeNe laser (with wavelength A = 0.6328 p,m> into the input fiber of a 50:50 coupler. The result is that the light from the laser is split in two parts. The fiber used was a nonpolarization-preserving monomode fiber with 4 brn core. A fiber polarization controller, realized by a few bends of variable orientation in the fiber [ 121, is used to maximize the visibility of the projected fringes. The ends of the fibers are mounted on a precision translator, so that the separation of the fibers (a) may be controlled to
Fig. 2. Optical geometry of the system.
micron resolution to obtain the desired fringe separation on the object. In order to have a set of parallel fringes, the fiber ends have been placed at the focal point (.f) of a collimating lens. Therefore, the fringe spacing ( p(,), on a screen (orthogonal to the illuminating axis), is given, to a first-order approximation, by ,‘,,
=
If‘
(1)
u
plane onto the photosensor of a CCD TV camera (we assume that the video signals from the TV camera are proportional to the light intensities). An interference filter, centered at 633 nm, is used to obviate the need of operating in dark ambient. The signal from the CCD camera is digitized in a frame grabber unit linked to a PC computer. If the vibration frequency is much higher than the framing frequency of the TV camera, the electrical video signal is proportional to the average intensity over one period,
so that, the resulting intensity profile can be shown to be tl II
(I( X,.v)), = 210/+ I + cos 2n;cps9 --P/W i [
where I,, is the intensity from each fiber end, and n is the refractive index of the medium (assumed homogeneous and equal to one). The optical geometry of the system is shown in Fig. 2. The phase term. present in Eq. (2). can be written in terms of the viewing axis. Indeed the X, coordinate can be written as (see Fig. 3) .r,=.~cost!-,-sint9.
I( X.!.)
I\
I fcos
This yields 2~UCOSB
(I(.r,J)),=K
I +.I,
i
[
hf
:,, tan29
I
(3)
If the surface under study is locally plane and is oscillating sinusoidally in the : direction ( : = ;,) CDSor) at any instant of time. the intensity in the object plane is
=2/,,
(5)
X(xfzacosortan-9)
2rr~.(.~-cos19+;,,coswfsin~9) I, L
Ii
(4)
The light intensity of the frrnge pattern can be converted to an electrical video signal by imaging the object
(6) where K is a constant. J,, is the zero-order Bessel function of the first kind, and ;o is the local amplitude of vibration. Eq. (6) represents the fringe pattern observed in real time on the video monitor. We note that the amplitude of the cosine term is modulated by a zero-order Bessel function whose argument is proportional to the local amplitude of vibration. The information about the local vibration amplitude can be demodulated by a Fourier transform method [13,141.
v ‘. I:..:1
L
!’
\.
Fig. 3. Geometry to write source axis coordinate in terms of the viewing axis.
20
G. Schirripa Spagnolo er al. /Optics
Communications
139 t 19971 I7-23
Fig. 4. Time average CCD images for the cantilever beam.
Now Eq. (6) can be rewritten as: (I(x,~)),=K,{1+5,(2~u, +J,(2n
distinct regions, corresponding to the three terms in Eq. (71, separated by a frequency uO,
zOtan6)exp(i2ru,x) uO zO tan*)
exp( - i2r
uO x)} ,
FT{(l(x,y)),}
=A(wJ)
+ C(u - uo>u)
(7) +C(--u-uoJ).
where
(9)
a cos6
u”=T’ The Fourier transform (FI’) of (I( x,y)),
(8) consists of three
Let us select with a filter mask the frequency region corresponding to C( u - uO,u) in Eq. (9). Next, we shift by an amount ug the selected region of
Fig. 5. TV monitor images of the /JoI functions, relative to Fig. 4. after FFT analysis. representing the maximum of I.I,,I and black representing the zero.
The images are codified
in gray levels, with white
G. Schirripa Spagmlo et al. / Optics Communications I39 (1997) 17-23
21
the spectrum in order to restore the function C(u,u). The inverse Fourier transform, with respect to u and v, yields B( x.y)
= K, J&77
u. z. tad),
(10)
where K, is a constant of proportionality. B(x, y) can be used to determine the vibration amplitude. It can be also shown on a TV monitor
I B( x,y)l a IJ,,(27r u. z,, tamY)l, so that the FFT provides
(11)
a high-quality
display
of lJ,(
fringes. To determine the constant K,, a FFT with the object in static position is realized. Generally this constant is function of the x, y coordinates. The local value of K, allows to minimize the effect of fringe and object surface curvature.
To confirm the usefulness of this method, some experiments were performed with the configuration shown
,
B
/’
---c /’
E”.a-
,/’ .’
$o.,
,I
’
/’
II
,,’ ,’
i
A/_.-
'E 0.4
O.Z-
/
,’
I
Fig. 7. Plot of the vibration amplitude of the cantilever a typical line from the fixed to the free end.
3. Experimental resdts
-A ....
I
,*;>
plate along
schematically in Fig. 1. First, a steel cantilever plate of dimension 50 mm by 40 mm by 0.25 mm thickness. was used as a test object. It was rigidly clamped at one end and free at the other. It was excited near the fundamental mode by acoustic coupling to a loudspeaker, driven by an audio oscillator. Fig. 4 shows some time-average CCD images. The images consist of the basic high-frequency “carrier fringe” pattern, modulated by bands of decreased visibility. These bands are centered on regions where the amplitude is such that the Bessel function in Eq. (11) is a minimum. The images shown in Fig. 4(A) , 4(B) and 4(C) were produced in the same way, except for the power of excitation changed to increase the amplitude of oscillation.
!
Fig. 6. Plot of the 1Jo1function. and vibration amplitude, averaged over the width of the cantilever plate, relative to Fig. 5. Amplitude at the free end of the cantilever is: (A) 0.5 mm; (Bl 0.8 mm; (0 1.1 mm.
Fig. 8. Time average CCD images for the plate
G. Schirripa Spagnoio et al. / Optics Communications
Fig. 9. Plot of the function 1J,,I obtained by m
Fig. 5 shows the images relative to Fig. 4 after FFT analysis. Fig. 6 is the plot of the IJ,I, and vibration amplitude, averaged over the width of the cantilever plate, relative to Fig. 5. Fig. 7 shows the vibration amplitude of the cantilever plate along a typical line from the fixed to the free end.
As a second test object a square steel plate, clamped along its boundary, having dimension 60 X 60 X 0.25 mm3 was used. Fig. 8 shows time-average CCD images of the plate vibrating near the fundamental mode. Fig. 9 is the plot of the function lJol obtained from the Fig. 8. The experiment has been realized with a fiber separation a = 0.25 mm, fringe separation on the target 2 mm, and viewing angle 8 = 60”. With our system we can evaluate changes in the J, function of about 5%, therefore the resolution achievable in the vibration amplitude measurement is about 50 km.
4. Conclusion A fiber optic projection interferometer and its application to vibration studies have been presented. This tech-
139 (I 997) I7-23
analysis of Fig. 8
nique combines the advantages normally associated with time-average holography with the convenience, low mechanical stability requirements, and real-time viewing of electronic speckle techniques, yet it produces a map of vibrational amplitude (not just nodal regions). Its sensitivity can be varied by adjusting the distance between the fiber ends, but in practice it is lower than the time average holography and works best for relatively large amplitude measurements. Indeed, the sensitivity can be adapted for vibration amplitudes in the range of 100 km to 10 mm. The electronic image system, the image digital processing and the fiber optics utilization eliminate many of the limitations inherent to the method proposed in Ref. [lo].
References [I] C.M. Vest, Holographic interferometry (Wiley, New York, 1979). [2] R.J. Pryputniewicz, Opt. Eng. 24 (1985) 843. [3] E. Vikhagen, Optics Comm. 69 (1989) 214. [41 C.S. Vikram, Study of vibration, in: Laser speckle and related phenomena. ed. P.K. Rastogi (Springer, Berlin, 1994) p. 293. [51 G. Schirripa Spagnolo, D. Paoletti, P. Zanetta. Optics Comm. 123 (1996) 41.
G. Schirripa Spagnolo et al. / Optics Communications [6] J.D. Valera, J.D.C. Jones. Electron. Lett. 29 (1993) 1789. [71 D.R. Burton, M.J. Lalor, Appl. Optics 33 (1994) 2939. [81 C. Quan, C.J. Tay. H.M. Shang, P.J. Bryanston-Cross. Optics Comm. I 19 (1995) 479. [91 D.P. Towers. C.H. Buckberry, B.C. Stockley. M.P. Jones. Meas. Sci. Technol. 6 (1995) 1242.
I39 (1997) 17-23
23
[lOI CM. Vest, D.W. Sweeney, Appt. Optics 1 I (1972) 449. I1 11 E. Hecht, Optics (Addison-Wesley, Reading, MA, 1987). [121 H.C. Lefevre, Electron Lett. 16 (1980) 778. [131 M. Takeda, H. Ina, S. Kobayashi, Opt. Eng. 31 (1992) 533. [I41 T. Kreis, J. Opt. Sot. Am. A 3 (1986) 847.
Optics Communications I39 (1997) 17-23
Vibration monitoring by fiber optic fringe projection and Fourier transform analysis Giuseppe Schirripa Spagnolo, Domenica Paoletti, Dario Ambrosini llipurtimento
di Energetica.
Uniwrsifir
drgli
Studi di L’Aquila.
Locali&
Mont&co
di Rob,
67040 Roio Poggio
- L’Aquila,
Italy
Received 14 November 1996; revised 15 January 1997; accepted 16 January 1997
Abstract A fiber optic fringe projection interferometer is proposed to obtain local amplitude vibration of a diffuse surface. A sinusoidal fringe pattern is projected onto the vibrating surface by an optical fiber version of Young’s interferometer. An image of this surface is captured by a CCD camera and fed to a PC through a frame grabber and hence displayed onto a TV monitor. The visibility of the interference pattern is modulated by a function of the local amplitude of vibration. Thus loci of constant vibrational amplitude can be directly observed as regions of low fringe visibility in the image of the surface. A quantitative analysis is performed by a Fast Fourier Transform (FFT) algorithm. Experimental results in agreement with the theory are presented. Keywords:
Vibration analysis; Fiber optic; Digital image analysis: Interferometry; Non-destructive testing
1. Introduction Many coherent optical techniques, based on holographic and speckle interferometry and its variants, have been devised for experimentally studying the vibrating surfaces [l-5]. Recently the introduction of digital image processing methods has led to new possibilities in this field, especially in the application of the speckle interferometry, but, despite any technical improvements, the limits to the quality of speckle images and the complex technologies of the holographic systems (the stringent stability requirements, etc.) are significant impediments to realize systems suitable for routine measurements in industrial environments. In this paper we propose to study vibrating objects by a sinusoidal fringe projection interferometer just developed for contour maps [6-91. The system is based on a periodic fringe pattern projected onto a vibrating surface which is imaged by a zoom lens onto the photosensor of a CCD camera. In most cases, the vibration frequency is much higher than the framing frequency of the CCD camera, so that the video signal, displayed on a monitor. is a time-averaged intensity pat-
tern. The resulting fringe pattern on the vibrating surface carries information regarding the local vibrational amplitude. Loci of constant vibrational amplitude can be directly observed as regions of high or low fringe visibility in the image of the surface. To get quantitative information on the local amplitude of vibration a Fast Fourier Transform (FIT) algorithm can be performed on the time-averaged CCD image. The FFI provides a high-quality display of J,, fringes that can be used to determine the vibration amplitude. The system is easy and cheap, with low mechanical stability requirements. It is well suited to measure vibrations with relatively large amplitude. The proposed method can be considered to some extent a new electronic version of the technique described by Vest and Sweeney [lo].
2. Experimental
setup and theory
The fringe projection system proposed in this paper is based on an optical fiber version of the well-known
0030-4018/97/$17.00 Copyright 0 1997 Elsevier Science B.V. All rights reserved. PII soo30-40 I X(97)00045-X
G. Schirripcr Spngnolo et al. / Optics Communications
Young’s interferometer [I I]. Interference fringes are formed with the coherent light of two single-mode fibers, which replace the pair of pin-holes in the Young’s configuration. In the classical version the power loss associated with the pin-holes is compensated by focusing the light onto their apertures. By doing so, however, Airy diffraction disks become visible and the regularity of the fringe pattern is altered as a consequence. These disadvantages are eliminated by the introduction of optical fibers, which allow Young’s fringes to be projected with a low power laser source.
139 11997) 17-23
The schematic diagram of the system is given in Fig. I. A graded-index (grin) lens is used to couple the light of a 3.5 mW HeNe laser (with wavelength A = 0.6328 p,m> into the input fiber of a 50:50 coupler. The result is that the light from the laser is split in two parts. The fiber used was a nonpolarization-preserving monomode fiber with 4 brn core. A fiber polarization controller, realized by a few bends of variable orientation in the fiber [ 121, is used to maximize the visibility of the projected fringes. The ends of the fibers are mounted on a precision translator, so that the separation of the fibers (a) may be controlled to
Fig. 2. Optical geometry of the system.
micron resolution to obtain the desired fringe separation on the object. In order to have a set of parallel fringes, the fiber ends have been placed at the focal point (.f) of a collimating lens. Therefore, the fringe spacing ( p(,), on a screen (orthogonal to the illuminating axis), is given, to a first-order approximation, by ,‘,,
=
If‘
(1)
u
plane onto the photosensor of a CCD TV camera (we assume that the video signals from the TV camera are proportional to the light intensities). An interference filter, centered at 633 nm, is used to obviate the need of operating in dark ambient. The signal from the CCD camera is digitized in a frame grabber unit linked to a PC computer. If the vibration frequency is much higher than the framing frequency of the TV camera, the electrical video signal is proportional to the average intensity over one period,
so that, the resulting intensity profile can be shown to be tl II
(I( X,.v)), = 210/+ I + cos 2n;cps9 --P/W i [
where I,, is the intensity from each fiber end, and n is the refractive index of the medium (assumed homogeneous and equal to one). The optical geometry of the system is shown in Fig. 2. The phase term. present in Eq. (2). can be written in terms of the viewing axis. Indeed the X, coordinate can be written as (see Fig. 3) .r,=.~cost!-,-sint9.
I( X.!.)
I\
I fcos
This yields 2~UCOSB
(I(.r,J)),=K
I +.I,
i
[
hf
:,, tan29
I
(3)
If the surface under study is locally plane and is oscillating sinusoidally in the : direction ( : = ;,) CDSor) at any instant of time. the intensity in the object plane is
=2/,,
(5)
X(xfzacosortan-9)
2rr~.(.~-cos19+;,,coswfsin~9) I, L
Ii
(4)
The light intensity of the frrnge pattern can be converted to an electrical video signal by imaging the object
(6) where K is a constant. J,, is the zero-order Bessel function of the first kind, and ;o is the local amplitude of vibration. Eq. (6) represents the fringe pattern observed in real time on the video monitor. We note that the amplitude of the cosine term is modulated by a zero-order Bessel function whose argument is proportional to the local amplitude of vibration. The information about the local vibration amplitude can be demodulated by a Fourier transform method [13,141.
v ‘. I:..:1
L
!’
\.
Fig. 3. Geometry to write source axis coordinate in terms of the viewing axis.
20
G. Schirripa Spagnolo er al. /Optics
Communications
139 t 19971 I7-23
Fig. 4. Time average CCD images for the cantilever beam.
Now Eq. (6) can be rewritten as: (I(x,~)),=K,{1+5,(2~u, +J,(2n
distinct regions, corresponding to the three terms in Eq. (71, separated by a frequency uO,
zOtan6)exp(i2ru,x) uO zO tan*)
exp( - i2r
uO x)} ,
FT{(l(x,y)),}
=A(wJ)
+ C(u - uo>u)
(7) +C(--u-uoJ).
where
(9)
a cos6
u”=T’ The Fourier transform (FI’) of (I( x,y)),
(8) consists of three
Let us select with a filter mask the frequency region corresponding to C( u - uO,u) in Eq. (9). Next, we shift by an amount ug the selected region of
Fig. 5. TV monitor images of the /JoI functions, relative to Fig. 4. after FFT analysis. representing the maximum of I.I,,I and black representing the zero.
The images are codified
in gray levels, with white
G. Schirripa Spagmlo et al. / Optics Communications I39 (1997) 17-23
21
the spectrum in order to restore the function C(u,u). The inverse Fourier transform, with respect to u and v, yields B( x.y)
= K, J&77
u. z. tad),
(10)
where K, is a constant of proportionality. B(x, y) can be used to determine the vibration amplitude. It can be also shown on a TV monitor
I B( x,y)l a IJ,,(27r u. z,, tamY)l, so that the FFT provides
(11)
a high-quality
display
of lJ,(
fringes. To determine the constant K,, a FFT with the object in static position is realized. Generally this constant is function of the x, y coordinates. The local value of K, allows to minimize the effect of fringe and object surface curvature.
To confirm the usefulness of this method, some experiments were performed with the configuration shown
,
B
/’
---c /’
E”.a-
,/’ .’
$o.,
,I
’
/’
II
,,’ ,’
i
A/_.-
'E 0.4
O.Z-
/
,’
I
Fig. 7. Plot of the vibration amplitude of the cantilever a typical line from the fixed to the free end.
3. Experimental resdts
-A ....
I
,*;>
plate along
schematically in Fig. 1. First, a steel cantilever plate of dimension 50 mm by 40 mm by 0.25 mm thickness. was used as a test object. It was rigidly clamped at one end and free at the other. It was excited near the fundamental mode by acoustic coupling to a loudspeaker, driven by an audio oscillator. Fig. 4 shows some time-average CCD images. The images consist of the basic high-frequency “carrier fringe” pattern, modulated by bands of decreased visibility. These bands are centered on regions where the amplitude is such that the Bessel function in Eq. (11) is a minimum. The images shown in Fig. 4(A) , 4(B) and 4(C) were produced in the same way, except for the power of excitation changed to increase the amplitude of oscillation.
!
Fig. 6. Plot of the 1Jo1function. and vibration amplitude, averaged over the width of the cantilever plate, relative to Fig. 5. Amplitude at the free end of the cantilever is: (A) 0.5 mm; (Bl 0.8 mm; (0 1.1 mm.
Fig. 8. Time average CCD images for the plate
G. Schirripa Spagnoio et al. / Optics Communications
Fig. 9. Plot of the function 1J,,I obtained by m
Fig. 5 shows the images relative to Fig. 4 after FFT analysis. Fig. 6 is the plot of the IJ,I, and vibration amplitude, averaged over the width of the cantilever plate, relative to Fig. 5. Fig. 7 shows the vibration amplitude of the cantilever plate along a typical line from the fixed to the free end.
As a second test object a square steel plate, clamped along its boundary, having dimension 60 X 60 X 0.25 mm3 was used. Fig. 8 shows time-average CCD images of the plate vibrating near the fundamental mode. Fig. 9 is the plot of the function lJol obtained from the Fig. 8. The experiment has been realized with a fiber separation a = 0.25 mm, fringe separation on the target 2 mm, and viewing angle 8 = 60”. With our system we can evaluate changes in the J, function of about 5%, therefore the resolution achievable in the vibration amplitude measurement is about 50 km.
4. Conclusion A fiber optic projection interferometer and its application to vibration studies have been presented. This tech-
139 (I 997) I7-23
analysis of Fig. 8
nique combines the advantages normally associated with time-average holography with the convenience, low mechanical stability requirements, and real-time viewing of electronic speckle techniques, yet it produces a map of vibrational amplitude (not just nodal regions). Its sensitivity can be varied by adjusting the distance between the fiber ends, but in practice it is lower than the time average holography and works best for relatively large amplitude measurements. Indeed, the sensitivity can be adapted for vibration amplitudes in the range of 100 km to 10 mm. The electronic image system, the image digital processing and the fiber optics utilization eliminate many of the limitations inherent to the method proposed in Ref. [lo].
References [I] C.M. Vest, Holographic interferometry (Wiley, New York, 1979). [2] R.J. Pryputniewicz, Opt. Eng. 24 (1985) 843. [3] E. Vikhagen, Optics Comm. 69 (1989) 214. [41 C.S. Vikram, Study of vibration, in: Laser speckle and related phenomena. ed. P.K. Rastogi (Springer, Berlin, 1994) p. 293. [51 G. Schirripa Spagnolo, D. Paoletti, P. Zanetta. Optics Comm. 123 (1996) 41.
G. Schirripa Spagnolo et al. / Optics Communications [6] J.D. Valera, J.D.C. Jones. Electron. Lett. 29 (1993) 1789. [71 D.R. Burton, M.J. Lalor, Appl. Optics 33 (1994) 2939. [81 C. Quan, C.J. Tay. H.M. Shang, P.J. Bryanston-Cross. Optics Comm. I 19 (1995) 479. [91 D.P. Towers. C.H. Buckberry, B.C. Stockley. M.P. Jones. Meas. Sci. Technol. 6 (1995) 1242.
I39 (1997) 17-23
23
[lOI CM. Vest, D.W. Sweeney, Appt. Optics 1 I (1972) 449. I1 11 E. Hecht, Optics (Addison-Wesley, Reading, MA, 1987). [121 H.C. Lefevre, Electron Lett. 16 (1980) 778. [131 M. Takeda, H. Ina, S. Kobayashi, Opt. Eng. 31 (1992) 533. [I41 T. Kreis, J. Opt. Sot. Am. A 3 (1986) 847.