Optical and Quantum Electronics (2005) 37:1199–1212 DOI 10.1007/s11082-005-4192-7

© Springer 2006

Dispersion compensation in optical coherence tomography with a prism in a rapid-scanning optical delay line meng-tsan tsai, i-jen hsu, chih- wei-lu, yih-ming wa n g , c h i a - w e i s u n , y e a n - w o e i k i a n g a n d c . c . y a n g∗ Department of Electrical Engineering, Graduate Institute of Electro-Optical Engineering, National Taiwan University, 1, Roosevelt Road, Section 4, Taipei, Taiwan, R.O.C. (∗ author for correspondence: E-Mail: [email protected]) Received 15 February 2005; accepted 15 September 2005 Abstract. We demonstrate the theoretical and experimental results of using a single prism in the rapid-scanning optical delay line of an optical coherence tomography (OCT) system for compensating the mismatches of the first- and second-order group delay dispersion (GDD) between the reference and sample arms. The analytical expressions for the first- and second-order GDD are derived based on the typically designed system configuration. Numerical results in varying various parameters are shown. An optimized set of parameters for efficient dispersion compensation in a practical fiber-based OCT system is obtained. The numerical result of the dispersion compensation is demonstrated. Also, the experimental implementation of such a dispersion compensation method is illustrated with the conditions similar to the numerical calculations. The compensation result is quite satisfactory. Key words: biophotonics, optical coherence tomograph, dispersion compensation, rapid-scanning optical delay line

1. Introduction Optical coherence tomography (OCT) has been proved a useful technique for noninvasive biomedical imaging (Huang et al. 1991). It can be used for in vivo real-time imaging with high spatial resolution. The developments of broadband super-continuum generation and extremely short pulses as the light sources of OCT systems have led to ultrahigh resolution around one micron in biological tissues (Kowalevicz et al. 2002; Marks et al. 2002; Povazay et al. 2002; Vabre et al. 2002 Bizheva et al. 2003; Unterhuber et al. 2003; Wang et al. 2003;). For high-speed scanning, a rapid-scanning optical delay line (RSODL) has been widely used in the reference arm of an OCT system (Rollins et al. 1998; Zvyagin et al. 2003). The typical configuration of an RSODL is based on the concept of Fourier-transform pulse shaping by placing a tilted mirror at the Fourier plane of a lens after a diffraction grating to produce a phase ramp in the frequency domain. This arrangement has the advantages of high-speed linear scan, group- and

1200

M.-T. TSAI ET AL.

phase-delay independencies, and feasible group-velocity dispersion compensation. In OCT implementations, fiber-based systems have been widely used because of the advantages of high flexibility and stability (Saxer et al. 2000; Guo et al. 2004; Todorovi´c et al. 2004). However, in using RSODL in the reference arm of a fiber-based OCT system, usually significant dispersion mismatch between the reference and sample arms is introduced, particularly when the fiber lengths of the two arms are different. Although methods of slightly adjusting the configuration of an RSODL can provide certain group-velocity dispersion compensation, the effect is limited. In particular, they usually cannot compensate the higher-order dispersion mismatch that will still broaden the interference fringe envelope and introduce side lobes. Such methods include one for real-time dispersion compensation by statically tilting the grating in the RSODL (Smith et al. 2002). However, this method can only compensate the dispersion that is linearly proportional to the scanning position for achieving depth-dependent (also called timedependent) dispersion compensation (including the first and second-order group delay dispersion). Also, by properly translating the position of the grating relative to the lens along the optical axis, the time-independent firstorder group delay dispersion (GDD) mismatch, which originates from the dispersion in fiber, can be partly compensated (Niblack et al. 2003). Nevertheless, even with the availability of the aforementioned methods, the timeindependent higher-order GDD from fiber cannot be compensated. Such dispersion mismatch issue becomes more important as the light source spectrum becomes broader. Other dispersion compensation methods have been reported. Dispersion correction with numerical algorithms has been widely implemented. A retrieval algorithm was proposed to reduce the side-lobes of the interference pattern such that the resolution could be effectively improved (Hsu et al. 2003). Numerical methods have been reported for compensating the depth-dependent sample dispersion (de Boer et al. 2001; Fercher et al. 2001, 2002; Marks et al. 2003a,b). Software dispersion compensation is quite useful in improving the axial resolution of an OCT system; however, it usually requires extra time for processing the image data and may slow down the imaging speed. Dispersion compensation was also realized by carefully controlling the fiber lengths in a fiber-based OCT system such that the dispersions in the reference and sample arms are exactly balanced (Chen and Li 2004). However, this method requires the precise control of fiber length down to the mm range and hence makes the system operation inflexible. Besides the dispersion compensation in the time-domain OCT, dispersion correction in the spectral domain OCT caught much attention (Wojtkowski et al. 2004). Recently, it was reported that dispersion compensation could be implemented with a prism pair or a piece of glass plate in

DISPERSION COMPENSATION IN OCT

1201

the reference arm (Kowalevicz et al. 2002; Bourquin et al. 2003; Herz et al. 2004). Although this method requires extra optical components and makes the system more complicated, it is quite effective in the OCT systems of extremely large spectral widths. In this paper, we demonstrate the detailed study on a simple method for depth-independent (time-independent) compensation of higher-order GDD mismatch in an OCT system. In this method, a prism (instead of a prism pair by Bourquin et al. (2003), or a piece of glass plate by Herz et al. (2004)) is placed between the grating and the lens of the RSODL. Because different optical paths are traced in different wavelengths within the prism and the refractive index of glass is dependent on wavelength in a nonlinear manner, compensation of the higher-order dispersion mismatch can thus be achieved. By properly designing the parameters of the prism, including the position, apex angle, orientation, etc., all the time-independent GDD (the first- and second-order GDD) can be compensated. Both theoretical and experimental results are demonstrated in this paper. This paper is organized as follows: In Section 2, analytical expressions for the first- and second-order GDD, based on ray tracing, are formulated. With the formulations, numerical results are discussed in Section 3. Then, experimental procedures and results are presented in Section 4. Finally, conclusions are drawn in Section 5.

2. Theories Figure 1 shows the configuration of an RSODL with a prism placed between the grating and the lens for dispersion compensation. The grating surface is arranged to be perpendicular to the diffracted signal of the center wavelength. The lens is placed such that the signals of different wavelengths can be parallel to each other right to the lens, and the signal of the center wavelength propagates along the optical axis of the lens toward the tilted mirror with its axis of rotation offset by x0 . The distance between the lens and the tilted mirror is equal to the focal length, f, of the lens. The wavelength-dependent phase delay can be calculated with ray tracing analysis and be expressed as (see Fig. 1) ψ(k, θ ) = k([OC] + [CDE] + [EG])

(1)

with [OC] = OA + nk AB + BC,

(2)

[EG] = nk EF + F G.

(3)

1202

M.-T. TSAI ET AL.

Fig. 1. The configuration of an RSODL with a prism inserted between the grating and the lens for dispersion compensation.

Here, [OC] denotes the optical path length between O to C, and OA denotes the physical distance between O and A. The notations nk and k stand for the refractive index of the prism and the wave number, respectively. The first- and second-order GDD are defined, respectively, as
d2 ψ

ψ  (k0 ) = (4) Dω ≡
dω2
c2 ω=ω0


d ψ

ψ  (k0 ) (1) = Dω ≡ c3 dω3
ω=ω0 3

(5)

Here,  is the phase delay, ω is the angular frequency, k0 is the wave number at the central wavelength, and c is the speed of light. Both GDD terms can be obtained by differentiating the wavelength-dependent phase delay with respect to the wave number k. With ray tracing analysis, the phase shift can be expressed as (k)] ψ(k) = 2k(Lr + [ODG] + H Q − Ls ) = 2k[δ + ψ

(6)

Here, δ ≡ Lr − Ls + H Q + [O  MO  ] ˜ ψ(k) = 2d1 cos β ⎡

  ⎤  2 X1 1 − X2 − X2 1 − X12 (n2 − sin2 φ)(1 − X2 ) 0 2 ⎦ +2d2 ⎣sin β + cos φ n20 − X22

(7)

DISPERSION COMPENSATION IN OCT



 2 2 +2(d3 − f ) X1 X2 + (1 − X1 )(1 − X2 )

  +2f θ (X1 1 − X22 − X2 1 − X12 ) − 2x0 θ  X1 ≡ sin γ n2k − sin2 (β + φ) − cos γ sin(β + φ)  X2 ≡ sin γ n20 − sin2 φ − cos γ sin φ

1203

(8) (9) (10)

After certain calculations, the first- and second-order GDD can be obtained, respectively, as Dω =

ψ  (k0 ) 4  2k0 = 2 ψ˜ (k0 ) + 2 ψ˜  (k0 ), 2 c c c

(11)

Dω(1) =

ψ  (k0 ) 6  2k0 = 3 ψ˜ (k0 ) + 3 ψ˜  (k0 ). c3 c c

(12)

Here, Lr and Ls are the path lengths of the reference (excluding RSODL) and sample arms. They are unimportant in our derivations. The notation d1 is the distance between the grating and the tip of the prism, d2 is the distance between the incident position of the light on the grating and the projection point of the prism tip on the grating, and d3 is the distance between the prism tip and the lens. Also, α and β are the incident and diffracted angles, respectively, of the input light on the grating. γ is the apex angle of the prism, φ is the angle between the grating and the prism, and θ is the tilt angle of the rotating mirror. Finally, n0 = 1 is the free-space refractive index. As mentioned earlier, the methods of translating and tilting the grating for dispersion compensation have been reported. Figure 2 shows the configuration of an RSODL without the prism. Again, the dispersion characteristics can be calculated with a ray tracing analysis. In this situation, the phase delay becomes ψ(k) = 2kδ − 4kx0 θ − 4kf θ sin β + 4k z cos β.

(13)

The first and second terms represent the dispersion effects produced in the RSODL. The third term results from the tilt of the grating and the fourth term originates from the translation of the grating. By differentiating the phase delay and using the grating equation, the first- and second-order GDD can be obtained. The contributions to the first- and second-order GDD, respectively, from the fourth term on the right-hand side of (13) are

1204

M.-T. TSAI ET AL.

Fig. 2. The configuration of an RSODL for dispersion compensation through tilting and translating the grating.

16π 2 m2 z c2 2 k03 cos2 θg  2mπ sin θg 48π 2 m2 z (1) 1+ Dω, z = 3 2 4 k0 cos2 θg c k0 cos2 θg

Dω, z = −

(14)

(15)

Here, θg is the tilt angle of the grating and z is the displacement of the grating from the focal plane of the lens. Also, m = −1 represents the diffraction order of the grating and is the grating period. If the dispersion is independent of the scanning angle θ, it is defined as time-independent dispersion. Since the dispersion corresponding to the first term of (13) is independent of θ, the effect of translating the grating, given in (14), can compensate the time-independent first-order dispersion. However, as shown in (15), it causes time-independent second-order GDD. Meanwhile, it cannot compensate time-dependent dispersion. The contributions to the first- and second-order GDD, respectively, from the third term on the right-hand side of (13) are 16π 2 m2 f θ sin θg c2 p2 k03 cos3 θg  2mπ sin θg 48π 2 m2 f θ sin θg (1) 1+ Dω,θ = 3 2 4 pk0 cos2 θg c p k0 cos3 θg Dω,θ = −

(16)

(17)

Here, one can see that although tilting the grating can introduce time-dependent dispersion for compensating the second term of (13) (also the dispersion from sample material), it generates time-dependent second-order GDD. The above discussions show that simply translating and tilting the

DISPERSION COMPENSATION IN OCT

1205

grating are insufficient for effective dispersion compensation. On the other hand, by using the configuration with the prism, as shown in Fig. 1, we can adjust those parameters, including d1 , d2 , d3 , φ, γ , etc., for compensating any time-independent, first- and second-order GDD in an OCT system, generated from the RSODL and the fiber length difference. Also, if the sample material dispersion can be estimated, those parameters can be optimized for efficiently compensating such a time-dependent dispersion mismatch. In the following, we will only discuss the details of time-independent dispersion compensation.

3. Numerical results In this section, we first numerically demonstrate the dispersion behaviors of the ROSDL with the prism by varying several key parameters. Such information will be used to optimize the parameters for compensating the dispersion in a practical OCT system. In Figs. 3–6, we show the numerical results of the first- and second-order GDD in the RSODL with the dispersion compensation prism by varying various parameters. In Fig. 3, the distance between the prism tip and the grating, d1 , is varied with other parameters fixed at d2 = 1.306 cm, d3 = 7.3 cm, f = 7.5 cm, n0 = 1, θ = 0.5◦ , α = 52.35◦ , = 2.5µm (400 lines/mm), γ = 45◦ , x0 = 1 mm, and φ = 10◦ . Here, one can see the slight increasing and decreasing trends with increasing wavelength for the first- and second-order GDD, respectively. Either order of GDD increases with increasing d1 value. Also, the considered d1 values mainly lead to normal dispersion. Figure 4 shows the similar data by varying the distance between the prism tip and the lens, d3 . Except that the d1 value is fixed at 1 cm and d3 is varied, all other parameters are the same as those leading to Fig. 3. Here, anomalous dispersion is obtained when d3 is smaller than 9 cm. Either the first- or second-order GDD increases with increasing d3 value. Figure 5 shows the similar results by varying the prism apex angle, γ . Here, except that d1 is fixed at 1 cm and γ is varied, all other parameters are the same as those for Fig. 3. One can see the broad distribution of either the first- or second-order GDD within the concerned wavelength range. Local maximum values of the first-order GDD exist for the cases of γ = 45◦ and 60◦ . Figure 6 shows the similar data by varying the grating period, . Here, the first- and second-order GDD have the similar wavelength-dependent behaviors. The dependence becomes stronger with increasing wavelength. For this figure, except that d1 is fixed at 1 cm and is varied, all other parameters are the same as those for Fig. 3. In Fig. 7, we compare the first- and second-order GDD between the RSODL and a 20 cm optical fiber. The dispersion values of the RSODL (dashed curves) are obtained by setting the parameters the same as those for Fig. 3 and d1 = 1 cm. The RSODL dispersion results are almost

1206

M.-T. TSAI ET AL.

Group delay dispersion (ps/nm)

0.25

(a)

d1=9 cm

0.20

7 cm 0.15

5 cm 0.10

3 cm 0.05

1 cm

0.00 600

700

800

900

1000

Second-order group delay dispersion 2 (fs/nm )

Wavelength (nm) 0.20

(b)

0.15 0.10

d1= 9 cm 7 cm 5 cm 3 cm 1 cm

0.05 0.00

600

700

800

900

1000

Wavelength (nm) Fig. 3. The first-order (part (a)) and the second-order (part (b)) GDD of the RSODL system with different d1 values. All other parameters are fixed.

independent of θ when its absolute value is smaller than 5◦ (the rotation limit in real operation). The fiber data (continuous curves) are based on those of a Corning fiber (model HI-780). The RSODL parameters were actually optimized to match the first- and second-order GDD simultaneously at 800 nm between the RSODL and the fiber. Therefore, if an OCT system consists of the designated RSODL (including the prism) and a 20 cm fiber length difference between the reference and sample arms, the dispersion mismatch in the whole OCT system can be effectively compensated around 800 nm. Figure 8 shows the calculation results of the OCT interference fringe envelopes in the cases of dispersion compensation with the prism (continuous curve) and the theoretical limit of zero dispersion mismatch (dashed curve). The results were obtained based on the assumption of a Gaussian

1207

DISPERSION COMPENSATION IN OCT

Group delay dispersion (ps/nm)

0.05

d3= 9 cm 7 cm

0.00

5 cm -0.05

3 cm 1 cm

-0.10

(a) -0.15 600

700

800

900

1000

Second-order group delay dispersion 2 (fs/nm )

Wavelength (nm)

0.1

d3= 9 cm 7 cm

0.0

5 cm

-0.1

3 cm -0.2

1 cm (b) -0.3 600

700

800

900

1000

Wavelength (nm) Fig. 4. The first-order (part (a)) and the second-order (part (b)) GDD of the RSODL system with different d3 values. All other parameters are fixed.

spectral shape with the central wavelength at 800 nm and the spectral FWHM at 60 nm. The OCT system consists of the aforementioned RSODL (with the prism) and the 20 cm fiber length difference. The dispersion of the RSODL was calculated by choosing θ = 0.5◦ . The result is insensitive to this choice. Without the dispersion compensation of the prism, the envelope FWHM is 35 µm (not shown). Due to the second-order GDD, the envelope is asymmetric. Through the dispersion compensation with the conditions shown in Fig. 7, the FWHM of the interference fringe envelope or the OCT resolution has been reduced to 5.17 µm, which is close to the theoretical limit of 4.71 µm. Because the second-order GDD is not completely compensated, the continuous curve in Fig. 8 is still slightly asymmetric. Also, some minor side lobes exist on the left-hand side of the envelope.

1208

M.-T. TSAI ET AL.

Group delay dispersion (ps/nm)

0.01

γ = 300

0.00

0

45

-0.01 -0.02 -0.03 -0.04

0

60

(a) -0.05 600

700

800

900

1000

Second-order grou delay dispersion 2 (fs/nm )

Wavelength (nm) 0.05

γ = 30o

0.00

o

45

-0.05 -0.10 -0.15

o

60

(b) -0.20 600

700

800

900

1000

Wavelength (nm) Fig. 5. The first-order (part (a)) and the second-order (part (b)) GDD of the RSODL system with different γ values. All other parameters are fixed.

4. Experimental procedures and results In the experiment, a fiber-based OCT system with an RSODL was built with a mode-locked Ti: Sapphire laser as the light source, as shown in Fig. 9. The laser spectral FWHM is 60 nm with the central wavelength at 800 nm. Therefore, the theoretical limit of the axial resolution is 4.71 µm. Corning fiber of model HI-780 is used. The fiber lengths in the sample and reference arms are 2 and 1.8 m, respectively. Therefore, the 20-cm fiber length difference causes dispersion mismatch besides that from the RSODL in the interfered signals. Figure 10 shows the interference fringe envelope obtained by scanning a glass plate surface with the grating in the RSODL appropriately translated and tilted. Here, one can see that the dispersion mismatch is far from effective compensation. In particular, the asymmetry shows the strong effect of the

1209

DISPERSION COMPENSATION IN OCT

(a)

Λ = 2.50 µ m

Group delay dispersion (ps/nm)

0.00

1.67 µ m -0.15

-0.30

1.25 µ m

-0.45 600

700

800

900

1000

Second-order group delay dispersion 2 (fs/nm )

Wavelength (nm) (b)

Λ = 2.50 µ m

0.0

1.67 µ m -1.5

-3.0

1.25 µ m

-4.5 600

700

800

900

1000

Wavelength (nm)

0.000

0.4

-0.002

0.3 0.2

-0.004

0.1 -0.006 0.0 -0.008 SMF (20cm) RSODL

-0.1

-0.010 600

700

800

900

1000

Second-order group delay dispersion 2 (fs/nm )

Group delay dispersion (ps/nm)

Fig. 6. The first-order (part (a)) and the second-order (part (b)) GDD of the RSODL system with different values. All other parameters are fixed.

Wavelength (nm)

Fig. 7. The first- and the second-order GDD of a typical piece of fiber of 20 cm in length (continuous curves) and the RSODL with a set of optimized parameters (dashed curves).

1210

M.-T. TSAI ET AL.

1.0

Normalized intensity

0.8 0.6 0.4

5.17 µm

4.71 µm

0.2 0.0 Dispersion compensation with prism Complete dispersion compensation

-0.2 -0.4 -20

-15

-10

-5 0 5 Wavelength (nm)

10

15

20

Fig. 8. Comparison of the calculated interference fringe envelope between the cases of dispersion compensation with the prism (continuous curve) and the theoretical limit of zero dispersion mismatch (dashed curve). The spectral FWHM of a Gaussian-shaped light source is assumed to be 60 nm centered at 800 nm.

Fig. 9. Experimental setup of the fiber-based OCT system. GSM: galvanometer scanning mirror; BS: fiber beam splitter; DAQ: data acquisition board.

second-order GDD. When a prism was inserted into the RSODL between the grating and the lens, as shown in Fig. 1, after the optimized adjustments, the FWHM of the interference fringe envelope was significantly reduced (down to 6.87 µm), as shown in Fig. 11. This fringe envelope is somewhat similar to the continuous curve in Fig. 8. Because the practical conditions cannot be exactly the same as the theoretical predictions, the fringe envelope width in experiment is slightly larger than the results in Fig. 8. However, the experimental result does show the effectiveness of the first- and second-order GDD compensations by using the prism in the RSODL.

DISPERSION COMPENSATION IN OCT

1211

Fig. 10. Experimental interference fringe envelope with dispersion compensation through translating and tilting the grating.

Fig. 11. Experimental interference fringe envelope with dispersion compensation through inserting a prism into the RSODL.

5. Conclusions In summary, we have demonstrated with theory and experiment the effectiveness of compensating the mismatches of the first- and second-order GDD between the reference and sample arms of an OCT system by inserting a single prism into the used RSODL. The analytical expressions for the first- and second-order GDD have been derived based on the typically designed system configuration. Numerical results of varying

1212

M.-T. TSAI ET AL.

different parameters were shown. An optimized set of parameters for effective dispersion compensation in a practical fiber-based OCT system was obtained. The numerical results of the dispersion compensation were demonstrated. Also, the experimental implementation of such a dispersion compensation method was illustrated with the conditions similar to the numerical calculations. The compensation result was quite satisfactory.

Acknowledgement This research was supported by National Health Research Institute, The Republic of China, under the grant of NHRI-EX93-9220EI.

References Bizheva, K., B. Povaˇzay, B. Hermann, H. Sattmann and W. Drexler. Opt. Lett. 28 707, 2003. Bourquin, S., A.D. Aguirre, I. Hartl, P. Hsiung, T.H. Ko and J.G. Fujimoto. Opt. Express 11 3290, 2003. Chen Y. and X. Li. Optics Express 12 5968, 2004. de Boer, J.F., C.E. Saxer and J.S. Nelson. Appl. Opt. 40 5787, 2001. Fercher, A.F., C.K. Hitzenberger, M. Sticker, R. Zawadzki, B. Karamata and T. Lasser. Opt. Express 9 610, 2001. Fercher, C.K. Hitzenberger, M. Sticker, R. Zawadzki, B. Karamata and T. Lasser Opt. Comm. 204 67, 2002. Guo, S., J. Zhang, L. Wang, J.S. Nelson and Z. Chen. Opt. Lett. 29 2025, 2004. Herz, P.R., Y. Chen, A.D. Aguirre and J.G. Fujimoto. Opt. Express 12 3532, 2004. Hsu, I.J., C.W. Sun, C.W. Lu, C.C. Yang, C.P. Chiang and C.W. Lin. Appl. Opt. 42 227, 2003. Huang, D., E.A. Swanson, C.P. Lin, J.S. Schuman, W.G. Stinson, W. Chang, M.R. Hee, T. Flotte, K. Gregory, C.A. Puliafito and J.G. Fujimoto. Science 254 1178, 1991. Kowalevicz, A.M., T. Ko, I. Hartl and J.G. Fujimoto. Opt. Express 10 349, 2002. Marks, D.L., A.L. Oldenburg and J.J. Reynolds. Opt. Lett. 27 2010, 2002. Marks, D.L., A.L. Oldenburg, J.J. Reynolds and S.A. Boppart. Appl. Opt. 42 204, 2003a. Marks, D.L., A.L. Oldenburg, J.J. Reynolds and S.A. Boppart. Appl. Opt. 42 3038, 2003b. Niblack, W.K., J.O. Schenk, B. Liu and M.E. Brezinski. Appl. Opt. 42 4115, 2003. Povazay, B., K. Bizheva, A. Unterhuber, B. Hermann, H. Sattmann, A.F. Fercher and W. Drexler. Opt. Lett. 27 1800, 2002. Rollins, A.M., M.D. Kulkarni, S. Yazdanfar, R. Ung-arunyawee and J.A. Izatt, Opt. Express 3 219, 1998. Saxer, C.E., J.F. de Boer, B. H. Park, Y. Zhao, Z.P. Chen and J.S. Nelson. Opt. Lett. 25 1355, 2000. Smith, E.D. J., A.V. Zvyagin and D.D. Sampson. Opt. Lett. 27 1998, 2002. Todorovi´c, M., S. Jiao and L.V. Wang. Opt. Lett. 29 2402, 2004. Unterhuber, A., B. Povaˇzay, B. Hermann, H. Sattmann and W. Drexler. Opt. Lett. 28 905, 2003. Vabre, L., A. Dubois and A.C. Boccara. Opt. Lett. 27 530, 2002. Wang, Y., Y. Zhao, J.S. Nelson and Z. Chen. Opt. Lett. 28 182, 2003. Wojtkowski, M., V.J. Srinivasan, T.H. Ko, J.G. Fujimoto, A. Kowalczyk and J.S. Duker. Opt. Express 12 2404, 2004. Zvyagin, A.V., E.D.J. Smith and D.D. Sampson. J. Opt. Soc. Am. A 20 333, 2003.