## Abstract

We present a fast calculation of the electromagnetic field near the focus of an objective with a high numerical aperture (NA). Instead of direct integration, the vectorial Debye diffraction integral is evaluated with the fast Fourier transform for calculating the electromagnetic field in the entire focal region. We generalize this concept with the chirp z transform for obtaining a flexible sampling grid and an additional gain in computation speed. Under the conditions for the validity of the Debye integral representation, our method yields the amplitude, phase and polarization of the focus field for an arbitrary paraxial input field on the objective. We present two case studies by calculating the focus fields of a 40×1.20 NA water immersion objective for different amplitude distributions of the input field, and a 100×1.45 NA oil immersion objective containing evanescent field contributions for both linearly and radially polarized input fields.

©2006 Optical Society of America

## 1. Introduction

The plane wave spectrum(PWS) method is a well-known and efficient technique for calculating the propagation and diffraction of electromagnetic (EM) fields. Its efficiency lies in the ability to propagate EM fields from one plane to another using the fast Fourier transform (FFT).

In microscopy this concept is the essence of the Debye approximation and is often used for the calculation of the EM field [1, 2, 3] near the focus of high numerical aperture (NA) objectives. Török et al. considerably expanded this concept for studying the focal field distribution and its distortions in stratified media commonly encountered in optical microscopy [5]. For a general and historical review on diffraction theory the reader is referred to Stamnes [6].

However, for focal field calculations in microscopy, in particular for optical systems with high NA, this classical problem turns into a computational challenge due to the highly oscillatory behavior of the involved functions. In addition, polarization effects cannot be neglected rendering this calculation long and tedious. Recent techniques in microscopy and tomography such as the extended focus field [7], microscopy beyond the Abbe resolution limit and pointspread function engineering as advanced by S. Hell and his group [8], or rigorous ab initio calculations for fluorescence fluctuation spectroscopy [9] amplify the demand for fast focal field calculations.

In this paper we revisit the Debye approximation and propose a novel and flexible implementation of the Debye integral incorporating the effects of amplitude, phase and polarization in an overall manner. This new implementation is particularly suited for rapid numerical evaluation and requires substantially less effort for calculating the amplitude, phase and polarization of an EM field distribution generated by a high NA microscope objective.

The organization of this paper is as follows: Section 2 introduces the Debye approximation, i.e. the general framework and formulae used in the remainder of this work. Section 3 outlines the implementation based on the fast Fourier transform (FFT) and establishes the sampling and border conditions for obtaining accurate numerical results. Finally, section 4 presents selected examples, firstly the calculation of the EMfield for a 40×1.20NA water immersionmicroscope objective, and secondly, for a 100×1.45 NA oil immersion objective taking into account the evanescent field contribution.

## 2. The Debye diffraction integral as Fourier transform

This section establishes the basic formalism based on the Debye diffraction integral and the formulation of this integral as a Fourier transform. The basic optical layout and the respective coordinate systems are shown in Fig. 1. We assume that this optical setup, i.e. the imaging system obeys Abbe’s sine condition (as usually fulfilled for microscope objectives).

A coherent, monochromatic wave field parallel to the optical axis crosses the aperture stop A, propagates towards the principal plane ℙ_{1} and is transferred to the principal plane ℙ_{2}. At ℙ_{2}, the wave field is refracted and focused towards the focal point **F**
_{2}. The point **P** lies on the principal plane ℙ_{2} and illustrates the focusing of a ray at ℙ_{2} towards the focal point **F**
_{2}. The spherical surface ℙ_{2} is centered at **F**
_{2} and the deflection angle *θ* at the position **P** is given by

where *r* is the off-axis coordinate of the incident wave, *R* the aperture stop radius, *NA* the numerical aperture of the objective and *n*_{t}
the index of refraction behind the ℙ_{2} surface. In our setup, the aperture A is placed in the back focal plane, which results in a telecentric imaging system.

Instead of the principal planes, pupils are frequently used for modeling the wave propagation through the objective. However, diffraction at the aperture stop *inside the objective* is not obvious if the incident wave is transferred directly from the entrance pupil to the exit pupil. Within our representation, the wave propagation from the aperture plane A to the principal plane ℙ_{1} is easily calculated with the PWS method or in most cases based on classical Fourier optics principles.

The incident field *E⃗*_{i}
(*r*,*ϕ*) at ℙ_{1} is decomposed into a radial component (p-polarized) and a tangential component (s-polarized). The unit vectors for p- and s-polarization are

where *ϕ* is the azimuth angle around the *z*-axis. Upon transmission, the unit vector *e⃗*_{p}
is deflected by *θ* and becomes

Hence, the amplitude, phase and polarization of the transmitted field at ℙ_{2} is

where *t*_{p}
(*θ*,*ϕ*) and *t*_{s}
(*θ*,*ϕ*) are the transmission coefficients (viz pupil function, apodization) for p- and s-polarization, respectively. Accumulated phase distortions, i.e. aberrations at ℙ_{2}, as well as attenuations, i.e. amplitude factors, are integrated in the complex parameters *t*_{p}
and *t*_{s}
. As we assume the incident field to be paraxial, the axial component *E*_{iz}
is small against the lateral components *E*_{ix,y}
and can be neglected even if the incident phase is not constant. In the Debye approximation, the transmitted field *E⃗*_{t}
is the *plane wave spectrum* of the focus field *E⃗* near **F**
_{2}. Hence, the electric field *E⃗* at a point (*x*, *y*, *z*) is obtained by integrating the propagated plane waves, viz

$$=-\frac{\mathit{if}}{\lambda 0}\underset{0}{\overset{\Theta}{\int}}\mathrm{sin}\theta \underset{0}{\overset{2\pi}{\int}}{\overrightarrow{E}}_{t}(\theta ,\varphi ){e}^{i\left({k}_{z}z-{k}_{x}x-{k}_{y}y\right)}d\varphi d\theta .$$

The phase factor
${e}^{{\mathit{ik}}_{{z}^{z}}}$
accounts for the phase accumulation when propagating along the *z*-axis, whereas the term
${e}^{-i\left({k}_{x}x+{k}_{y}y\right)}$
represents the phase difference of the wave front at off-axis points (*x*, *y*, *z*) with respect to the on-axis point (0, 0, *z*). The integration extends over the solid angle Ω under which ℙ_{2} is observed at **F**
_{2}, i.e. sinΘ=*NA*/*n*_{t}
. The wave vector *k⃗*_{t}
is given in spherical coordinates *θ* and *ϕ* by

The evaluation of Eq. (5) is usually performed with a direct numerical integration taking into account the coordinate transformations, which results in the Richard-Wolf integral representation [2, 3]. Instead of the common ansatz, a (*θ*,*ϕ*)-sampling keeping dΩ=sin*θ*d*θ*d*ϕ* constant is obtained by using cos*θ*_{m}
=1-*m*ΔΘ with *m*∊ℕ. For *m*∊{1…*M*} and *n*∊{1…*N*}, the sampling grid is defined by

At *θ*=0, a sampling point with a weight of dΩ=${\pi \theta}_{1}^{2}$/4 is added. Besides minimizing the number of sampling points along *θ*, the calculation of the integrand and its integration can be merged in a single matrix product resulting in a further reduction of the computation time [4].

The outlined evaluation of the Debye diffraction integral (5) is quite fast, but still much slower than the conventional computation of a Fraunhofer diffraction integral. However, Eq. (5) can be easily rewritten as a Fourier transform by splitting the phase factor into a lateral and an axial term, and by performing the integration over ℙ_{1} instead of ℙ_{2}. Using Eq. (1) and (6), the integration step dΩ for a sampling over ℙ_{2} is projected onto ℙ_{1}, which yields

Insertion of this sampling step into Eq. (5) results in

Extending now the integration over (*k*_{x}
, *k*_{y}
)∊ℝ^{2} by setting |*E⃗*_{t}
|=0 for *r*>*R* allows to rewrite the Debye diffraction integral as a Fourier transform of the weighted field *E⃗*_{t}
, which finally results in

This is the main result of this work. The Debye integral is now expressed as a Fourier transform of the field distribution in the aperture A. The similarity of this expression with the conventional Fraunhofer diffraction integral is obvious. For a low NA imaging system, the weighting factor is approximated by 1/cos*θ*≈1 and Eq. (10) is equivalent to the Fraunhofer diffraction integral.

## 3. Numerical implementation

The numerical implementation is straightforward. A fast Fourier transform (FFT) of the weighted field at ℙ_{2} is used for the numerical evaluation of Eq. (10). For an equidistant sampling *k*_{x}
=*m*Δ*K* and *k*_{y}
=*n*Δ*K* with Δ*K*=*k*
_{0}
*NA*/*M*, viz *M* sampling points over the aperture radius, the sampling points on ℙ_{2} are

Multiplication of the integration step (Δ*K*)^{2} with the prefactor of Eq. (10) yields the numerical implementation of Eq. (10) as

Typically, the FFT is more than 100× faster than the direct integration of Eq. (5) with matrix multiplication. A good accuracy is achieved for 4*M*
^{2} ≳ 100×100 sampling points over Ω, but care has to be taken in order to avoid artifacts due to sampling and aliasing. Subsequently, the necessary conditions for obtaining accurate results are investigated [10].

#### 3.1. Sampling condition

The propagation factor
${e}^{{\mathit{ik}}_{{z}^{z}}}$
in Eq. (10) has to be calculated with high resolution for accurate results [11]. This imposes a condition on the phase discretization, i.e. the phase *k*_{z}*z* must not change by more than *π* between neighboring sampling points in the aperture plane A. With
${k}_{z}=\sqrt{{k}_{t}^{2}-{k}_{\mathit{xy}}^{2}}$
, the sampling condition can be expressed as

where Δ*K*=*k*
_{0}
*NA*/*M* and
$max\mid \mathrm{tan}\theta \mid =\frac{\mathit{NA}}{\sqrt{{n}_{t}^{2}-{\mathit{NA}}^{2}}.}$
. This immediately leads to a condition for the minimum number of sampling points

solely determined by the system parameters. For the numerical evaluation, an oversampling of about 3× is sufficient for improving the accuracy of the result. In addition, a lower limit of *M*≳50 reveals necessary for an accurate sampling of *ϕ*. Deviations from these sampling conditions result in granular artifacts as seen in Fig. 3(a). As a typical value for *M*, we have chosen *M*=125 for the focus field calculation of a 1.20 NA water immersion objective (see the example 4.1). A high accuracy is obtained for |*z*|≲25*λ*
_{0}, corresponding to ≈12 *µ*m at a wavelength of 488 nm.

#### 3.2. Sampling step

The focus field *E⃗* is obtained for the sampling positions (*m*Δ*x*,*n*Δ*y*, *z*). With Δ*k*=2*π*/*N*Δ*r* and Δ*r*=*f*Δ*K*/*k*_{t}
, the sampling step in the *xy*-plane is

where *N*>4*M* is the number of FFT sampling points per transformed dimension (see also Fig. 2, where the arrows span over 2*M*+1 samples and the padded dimension over *N* samples). For optimal FFT performance, it is best to set *N*=2
^{s}
with *s*∈ℕ. Respecting the condition (14), *M* can be adjusted to fit Δ*x* and Δ*y*. Along the *z*-direction, the sampling can be chosen arbitrarily by respecting the limits given above.

#### 3.3. Aliasing suppression

Due to the Debye diffraction integral expressed in Eq. (10), the field *E⃗*_{t}
is the plane wave spectrum of the focus field *E⃗*. Usually, the smallest area (aperture matrix) containing *E⃗*_{t}
≠0 is transformed (see Fig. 2). The spectral product
${e}^{{\mathit{ik}}_{{z}^{z}}}\times \frac{{\overrightarrow{E}}_{t}}{\mathrm{cos}\theta}$
in Eq. (10) represents a spatial convolution
$\overrightarrow{E}=\U0001d4d5\left({e}^{{\mathit{ik}}_{{z}^{z}}}\right)*\U0001d4d5\left(\frac{{\overrightarrow{E}}_{t}}{\mathrm{cos}\theta}\right)$
. In general, the result of the convolution is non-zero on an area larger than the aperture size, which may cause aliasing [12]. Therefore, the aperture matrix is enlarged by zero padding to at least twice its dimensions before performing the transform. In a final step, simple cropping of the transform output removes the padding.

Because we are only interested in the field near the focus, typically over a range of several wavelengths, we limit the computation of the FFT to this region of interest (Fig. 2). The transmitted field *E⃗*_{t}
is padded with zeros along the first dimension. In this dimension, the FFT is calculated and the result cropped. Along the second dimension, the same procedure is applied on the intermediate result. Zero padding simultaneously suppresses aliasing and refines the sampling grid for the focus field. Using two one-dimensional FFTs with intermediate cropping and zero padding minimizes the numerical processing cost.

#### 3.4. Aperture rim smoothing

Figure 3 shows the spectra log|FFT(*U*)| for a circular aperture with radius *R*. As already stated, the field *U* vanishes outside the aperture for *r*>*R*, whereas inside the aperture for *r*<*R*, the field is given as *U*=*U*
_{0}. This discretization leads to a serrated aperture rim inducing granular artifacts at higher frequencies. Hence, the expected Airy function is only seen at low frequencies (central region in Fig. 3(a), please note the logarithmic scale). A smooth sampling of the aperture rim improves the accuracy of the spectrum [13]. In Fig. (Fig. 3(b)) the rim was sampled with the hyperbolic tangent as

where Δ*R*=*R*/30 was the sampling step. The granular artifacts are efficiently reduced and the FFT approximates the Airy function with a good accuracy over a much larger area. Figure 4 shows a comparison of cross-sections of the spectra on the meridian *k*_{y}
=0. Overall, for values |*k*_{x}
|>60×2*π*/256Δ*R*, the ‘sharp’ spectrum shows granular artifacts, whereas the ‘smooth’ spectrum approximates well the Airy function.

#### 3.5. Generalization based on the chirp z transform

We demonstrated the importance of zero padding while respecting the sampling condition (14). These constraints led to a minimal number of sampling points *N*=2^{s} for the FFT (*s*∈ℕ). The corresponding number of sampling points *M* over the aperture radius often exceeds the initial guess based Eq. (14). In such cases, the chirp z transform (CZT) is computationally faster than the FFT. In summary, the CZT (a) allows breaking the relationship between *M* and *N*, (b) allows an implicit frequency offset, and (c) internalizes the zero padding. Applying this generalization, we adapted the sampling step in the focus field independently of the sampling step in the input field, introduced an additional shift of the region of interest, and finally improved the computational efficiency.

Let *z*_{m}
∀*m*∈[0,*M*-1] be a discrete representation of a spatial signal *z*(*r*=*m*Δ*r*). The discrete Fourier transform (DFT) at a frequency *k*=*n*Δ*k*∀*n*∈[0,*N*-1] is then obtained with

The FFT is a particular case of the DFT with Δ*k*=2*π*/*M*Δ*r* and *N*=*M*. For Δ*k*<2*π*/*M*Δ*r*, a zero padding is implicitly contained in Eq. (17). Comparing the DFT with the CZT defined by

yields *a*=1 and *w*=*e*
^{-iΔk} for obtaining the DFT as a particular case of the general CZT. Setting
$a={e}^{{\mathit{ik}}_{0}}$
shifts the frequency domain by *k*
_{0} (see above). Furthermore, Eq. (18) can be rewritten as a convolution

that can be evaluated using two (*M*+ *N*-1) point FFTs (a third one can be precomputed) [14].

Based on the CZT, our computationmethod can be extended for low NA systems or for focus fields with a large axial span. In such cases, the sampling grid becomes distorted over the focus depth [15, 16]. But within the framework of the CZT, this distortion can be compensated by a non-linear scaling proportional to the effective NA under which the aperture A is observed at ℙ_{2} from a point (0,0,*z*) on the axis. As a result, the sampling Δ*k* depends upon the axial position *z*, i.e. Δ*k*(*z*)=Δ*k*(0)*f*/(*f*+*z*) with Δ*k*(0)=Δ*k* as defined before. Using the CZT, the additional calculations remain restricted to the repeated computation of
$\mathrm{FFT}\left({w}^{\frac{-{m}^{2}}{2}}\right)$
because *w* varies now with *z*.

## 4. Selected examples

This section presents example calculations for two high NA microscope objectives. In the first example of a 1.20 NA water immersion objective, the variation of different amplitude distributions (apodization) in the aperture A are discussed. For the second example, a 1.45 NA oil immersion objective was chosen as used in total internal reflection microscopy. The refraction at a cover glass-water interface at the focus is added and the effect of different polarization distributions in the aperture plane A are discussed.

Before presenting these specific examples, the transmission coefficients *t*_{p}
and *t*_{s}
between the principal planes ℙ_{1} and ℙ_{2} need to be defined. We present the microscope objective as an optical system of only 2 optical interfaces and a convex interface into the immersion medium *n*_{t}
. To this end, the three interfaces provide a physical model for deflection angles *θ*∈[0,*π*/2). The amplitude transmission efficiency, i.e. apodization, and the polarization are obtained based on the Fresnel equations.

If the glass lens has an index of refraction *n*_{g}
and the immersion medium *n*_{t}
, the Fresnel transmission coefficients are calculated for the succession of the air(*n*_{a}
)-glass(*n*_{g}
)-air(*n*_{a}
)-immersion(*n*_{t}
) interfaces. The corresponding deflection angle *θ*_{ij}
at each interface was chosen proportional to the difference of the index of refraction, viz *θ*_{ij}
∝|*n*_{i}
-*n*_{j}
|. With *n*_{a}
=1, the Fresnel transmission coefficients are then

for p-polarization and

for s-polarization, respectively.

#### 4.1. 1.20 NA water immersion objective

Figure 5 shows the focus intensity for a nearly uniform and a Gaussian illumination in the back aperture of a 1.20 NA water immersion objective. For Δ*x*=Δ*y*=20 nm, Δ*z*=50 nm and *M*=100, a 2.0 GHz Pentium 4 processor computed the field within a volume of 3 *µ*m×3 *µ*m× 5 *µ*m i.e. 150 ×150 ×100 sampling points in less than 40 seconds. Taking the symmetry into account, the volume was further extended to 6 *µ*m×6 *µ*m×10 *µ*m.

In Fig. 5(a), the aperture was overfilled and the resulting focus field shows the well-known symmetry break of vectorial focus fields, for comparison the Airy profile was added. In Fig. 5(b), the aperture was underfilled to about 60% and the field becomes approximately gaussian. Figure 6 shows the electric fields along two major axes through the focus. For an overfilled aperture, the Airy profile (based on a scalar, paraxial approximation) is a good estimation of the electric field along the *y*-axis. For an underfilled aperture, the diameter of the central lobe is ≈25% larger but the side lobes vanish quickly. In both cases, the polarization leads to a larger x-extension compared to the *y*-extension.

Figure 7 and 8 show the intensity on the major planes through the focus. The polarization dependent extensions of the lobes along the major axes *x* and *y* creates a transition zone where the fringe contrast is diminished.

#### 4.2. 1.45 NA oil immersion objective

As a second example, we calculate the focus field of an objective designed for total internal reflection fluorescence (TIRF). The objective uses immersion oil with an index of refraction matching the cover slip. Its NA of 1.45 is higher than the index of refraction of the sample (*n*_{s}
=1.33, aqueous solution). This generates a partially evanescent focus field at the cover slip–sample interface. Depending upon the illumination of the aperture, the focus field can be fully propagating or fully evanescent. A fully propagating field can be calculated easily with the procedure outlined above. However, the evanescent field contribution needs an additional consideration for obtaining the total focus field.

First we determine the plane wave spectrum *E⃗*_{t}
at the immersion oil-cover slip interface. Next, the refraction at this interface and the cover slip–sample interface is calculated in order to obtain the plane wave spectrum *E⃗*_{s}
in the sample (water). Finally, applying the Fourier transform on the weighted and propagated spectrum
$\frac{{e}^{{\mathit{ik}}_{{z}^{z}}}{\overrightarrow{E}}_{s}}{\mathrm{cos}\theta}$
yields the focus field. As before, the angle *θ* and the weighting factor 1/cos*θ* are calculated in the immersion oil. But concerning the sampling condition, a specific issue related to the cover slip–sample interface (14) needs to be considered. The highest angles *θ* result in total internal reflection at the cover slip–sample interface. At the critical angle *θ*_{c}
=arcsin(*n*_{s}
/*n*_{t}
), *k*_{z}
vanishes. For higher angles, *k*_{z}
takes an imaginary value and the sampling condition (14) is relaxed because
${e}^{{\mathit{ik}}_{{z}^{z}}}$
becomes just an amplitude factor. The problem arises at *θ*_{c}
where the sampling condition (13) results in a singularity. Let *M*′ be the number of sampling points over *θ*<*θ*_{c}
. For avoiding this singularity at *θ*_{c}
, the sampling is chosen such that (*M*′+1/2)Δ*K*=*k*_{s}
, i.e. *θ*_{c}
falls between two sampling points. Inserting *M*=(*M*′+1/2)*NA*/*n*_{s}
,

into Eq. (13) then yields a generalized sampling condition

A 7× oversampling is used for improving the accuracy of the result, in particular at off-axis points. In addition, a lower limit of *M*≳100 was used for |*z*|→0.

Because the field is calculated in the sample space, *k⃗*_{t}
is replaced by

where *n*_{s}
sin*θ*′=*n*_{t}
sin*θ*. The unit vector *e⃗*_{r}
for p-polarization becomes

Figure 9 shows the focus field of a 100×1.45NA oil immersion objective. The aperture of the objective was overfilled, resulting in a partially evanescent field at the focus, where the cover slip–sample (water) interface was placed. As for the former example, the central lobe extends less in the *y*- than the *x*-direction for linear polarization (Fig. 9(a)). The focal volume is reduced to about 1/7 compared to the former water immersion objective. Selecting a radially polarized input field results in a rotationally symmetric focus field as shown in Fig. 9(b). On the optical axis, the electric field becomes purely *z*-polarized. For a distance *z*≲0.3 *µ*m, this *z*-component is dominant. Further away from the cover slip–sample interface, the *xy*-components prevail, which results in an annular field distribution.

The fine structure at the interface is due to the evanescent wave contribution with incidence angles above the critical angle. For instance, Fig. 12 shows the weighted field *E*_{s}
/cos*θ* for the linear polarization. At the critical angle (*NA*=1.33), the field amplitude triples in the *x*-direction and doubles in the y-direction, respectively, hence marking the abrupt transition from propagating to evanescent fields.

## 5. Conclusions

We showed a fast and simple implementation of the vectorial Debye integral for calculating the focus field of high NA objectives for arbitrary amplitude, phase and polarization distributions of the input field. The numerical evaluation with the fast Fourier transformis extremely fast and allows a high flexibility of the input field. The result is accurate under the conditions for the validity of the Debye integral representation of focused fields [17, 18] and the given sampling conditions. For low NA, it converges quite naturally to a focus field given by the Fraunhofer approximation. With the chirp z transform, we extended our calculations to low NA focus fields requesting a non-linear scaling as shown by Li and Hsu [15, 16]. Table 1 summarizes the performance of the different calculation methods on a personal computer.

In addition, we used a generalized pupil function (apodization) of high NA objectives taking into account amplitude and polarization distributions. The pupil function incorporates wave front aberrations as contained in real objectives as well as Fresnel transmission coefficients. Based on these Fresnel coefficients, it is straightforward to include wave propagation through stratified media.

In summary, our method allows fast and accurate calculations of the focus field in the entire focal region, which opens the path to fast simulations for point spread function engineering and image deconvolution in three-dimensional light microscopy.

## Acknowledgements

We are grateful to Herbert Gross, Carl Zeiss Oberkochen for many valuable comments and discussions. The support of the Swiss National Science Foundation (SNSF) (contract number 200021-103333) is greatly acknowledged.

## References and links

**1. **P. Debye, “Das Verhalten von Lichtwellen in der Nähe eines Brennpunktes oder einer Brennlinie,” Ann. Phys. **30**, 755–776 (1909). [CrossRef]

**2. **E. Wolf, “Electromagnetic diffraction in optical systems, I. An integral representation of the image field,” Proc. R. Soc. London Ser. A **253**, 349–357 (1959). [CrossRef]

**3. **B. Richards and E. Wolf, “Electromagnetic diffraction in optical systems, II. Structure of the image field in an aplanatic system,” Proc. R. Soc. London Ser. A **253**, 358–379 (1959). [CrossRef]

**4. **
Typically, a good accuracy is achieved for *M*≳50 and *N*≳200 sampling points.

**5. **P. Török and P. Varga, “Electromagnetic diffraction of light focused through a stratified medium,” Appl. Opt. **36**, 2305–2312 (1997). [CrossRef] [PubMed]

**6. **J.J. Stamnes, *Waves in Focal Regions: propagation, diffraction and focusing of light, sound and water waves*, Hilger, Bristol UK (1986).

**7. **G. Mikula, A. Kolodziejczyk, M. Makowski, C. Prokopowicz, and M. Sypek, “Diffractive elements for imaging with extended depth of focus,” Opt. Eng. **44**, 058001 (2005). [CrossRef]

**8. **N. Huse, A. Schönle, and S.W. Hell, “Z-polarized confocal microscopy,” J. Biomed. Opt. **6**, 480–484 (2001). [CrossRef]

**9. **J. Enderlein, I. Gregor, D. Patra, T. Dertinger, and U.B. Kaupp, “Performance of Fluorescence Correlation Spectroscopy for Measuring Diffusion and Concentration,” Chem. PhysChem. **6**, 2324–2336 (2005). [CrossRef]

**10. **
For simplification, the sample indices *kl* and *mn* will be omitted further on.

**11. **M. Mansuripur, “Certain computational aspects of vector diffraction problems,” J. Opt. Soc. Am. A **6**, 786–805 (1989). [CrossRef]

**12. **M. Sypek, “Light propagation in the Fresnel region. New numerical approach,” Opt. Comm. **116**, 43–48 (1995). [CrossRef]

**13. **P. Luchini, “Two-dimensional numerical integration using a square mesh,” Comp. Phys. Comm. **31**, 303–310 (1984). [CrossRef]

**14. **J. L. Bakx, “Efficient computation of optical disk readout by use of the chirp z transform,” Appl. Opt. **41**, 4897–4903 (2002). [CrossRef] [PubMed]

**15. **Y. Li and E. Wolf, “Three-dimensional intensity distribution near the focus in systems of different Fresnel numbers,” J. Opt. Soc. Am. A **1**, 801–808 (1984). [CrossRef]

**16. **W. Hsu and R. Barakat, “Stratton-Chu vectorial diffraction of electromagnetic fields by apertures with application to small-Fresnel-number systems,” J. Opt. Soc. Am. A **11**, 623–629 (1994). [CrossRef]

**17. **E. Wolf and Y. Li, “Conditions for the validity of the Debye integral representation of focused fields,” Opt. Comm. **39**, 205–210 (1981). [CrossRef]

**18. **P. Török, “Focusing of electromagnetic waves through a dielectric interface by lenses of finite Fresnel number,” J. Opt. Soc. Am. A **15**, 3009–3015 (1998). [CrossRef]