## Abstract

An algorithm is reported for the design of a phase-only diffractive optical element (DOE) that reshapes a beam focused using a high numerical aperture (*NA*) lens. The vector diffraction integrals are used to relate the field distributions in the DOE plane and focal plane. The integrals are evaluated using the chirp-z transform and computed iteratively within the Method of Generalized Projections (MGP) to identify a solution that simultaneously satisfies the beam shaping and DOE constraints. The algorithm is applied to design a DOE that transforms a circularly apodized flat-top beam of wavelength *λ* to a square irradiance pattern when focused using a 1.4-*NA* objective. A DOE profile is identified that generates a 50*λ*×50*λ* square irradiance pattern having 7% uniformity error and 74.5% diffraction efficiency (fraction of focused power). The diffraction efficiency and uniformity decrease as the size of the focused profile is reduced toward the diffraction limited spot size. These observations can be understood as a manifestation of the uncertainty principle.

©2008 Optical Society of America

## 1. Introduction

Beam shaping diffractive optical elements (DOEs) are becoming an essential part of many optical systems. A DOE can be used to modify the amplitude, phase, and polarization of an incident beam so that it focuses into a targeted irradiance distribution at the image plane [1–3]. Beam shaping can enhance performance in optical lithography, laser-based materials processing, direct laser writing, surgical applications, and optical data storage [4]. For many applications the optimum irradiance distribution consists of a flat-top profile having a defined geometry within the focal plane. Such irradiance patterns can be characterized in terms of the diffraction efficiency, *κ*, and uniformity error, *δ*. The diffraction efficiency quantifies the fraction of total optical power directed into the targeted region of interest and the uniformity error provides a measure of flatness in the irradiance distribution across that region.

Many excellent scalar techniques have been reported for designing beam shaping DOEs. These approaches are based on methods that include geometric mapping [5, 6], analytical solution [7], iterative processes [8–11], and genetic optimization [12]. Although exceptional results have been achieved with these algorithms, they are all based on scalar diffraction theory and as such are only valid in the paraxial domain of diffractive optics [13]. For systems with high numerical aperture (*NA*), depolarization effects are significant [14], so vectorial diffraction theory must be used in the DOE design process. This becomes particularly challenging because the overall beam shape is determined by the summed irradiance of the *x*-, *y*-, and *z*-polarized electric fields. Although the field components are orthogonal, they are not entirely independent because each is reshaped by a common DOE. As a result, the DOE must be designed to reshape the field components collectively so the irradiance of the *total* field matches the targeted beam shape. One report of vectorial beam shaping has appeared, but the method was not applied to high-NA systems [15]. Given that high-*NA* systems are being increasingly employed in frontier technologies, further applications of beam shaping will be stymied unless accurate methods for vectorial beam shaping are developed.

In this work we report a vectorial beam shaping algorithm that can be used to design phase-only DOEs for use under high-*NA* conditions. The algorithm was developed by incorporating the vector diffraction integrals [16] into the Method of Generalized Projections (MGP) [17]. The diffraction integrals are used to interrelate the DOE phase profile and the resulting vectorial electric field in the focal plane. The integrals are evaluated using the chirp-z transform [18] to improve computational speed and accuracy. Iterative projection of constraints in the pupil and focal planes progressively forces the simulation toward a DOE phase profile that generates the targeted beam shape. The new algorithm is applied to the problem of designing a phase-only DOE that transforms a circularly apodized flat-top input beam into a square flat-top irradiance distribution when focused using a 1.4-*NA* objective. In beam shaping, high diffraction efficiency and low uniformity error are known to be mutually exclusive characteristics that must be considered jointly in optimizing DOEs [19]. In this work, we also investigate how *κ* and *δ* change as the size of the focused beam profile approaches the diffraction-limited spot size.

## 2. Theory

#### 2.1. Vector diffraction integrals

The optical geometry of the focusing system is depicted in Fig. 1. A DOE and an aberration free aplanatic lens which fulfils the sine condition and has focal length *f* are positioned such that their optical axes are collinear with the *z*-axis of a Cartesian coordinate system whose origin is located at the Gaussian focus of the lens. The numerical aperture of the lens is *NA*=1.4 in all calculations. Monochromatic linearly polarized plane waves, with electric field vector parallel to the *x*-axis, propagate along the *z*-axis, passing through the DOE and entering the pupil of the lens. The light focuses into a medium of refractive index *n*=1.516. In the absence of the DOE, this situation is consistent with common applications of high-*NA* oilimmersion objective lenses.

The electric field at an arbitrary point *P*(*x _{f}*,

*y*) in the focal plane (

_{f}*z*=0) can be calculated using the vector diffraction integrals [16, 20] given as

_{f}At any point within the optical system the irradiance is *I*=(1/2)*ncε _{0}*|

*|*

**E**^{2}. The reshaped beam is the spatial map of the focused irradiance

*I*(

_{f}*x*,

_{f}*y*) for all

_{f}*P*. The speed of light and electric permittivity in vacuum are

*c*and

*ε*, respectively. The wave number of light transmitted through the lens is

_{0}*k*=2

_{t}*π*/

*λ*=[

*k*

^{2}

*+*

_{x}*k*

^{2}

*+*

_{y}*k*

^{2}

*]*

_{z}^{1/2}, with

*k*,

_{x}*k*, and

_{y}*k*being the plane wave components, and

_{z}*λ*is the wavelength within the medium. The

*NA*of the lens system sets

*k*

_{max}=

*k*/

_{t}NA*n*. The function

*T*(

*k*,

_{x}*k*)exp[

_{y}*iΦ*(

*k*,

_{x}*k*)] describes the transmission amplitude (

_{y}*T*) and phase (

*Φ*) of the DOE. The amplitude of the incident electric field

*E*is assumed to be spatially constant, so this term was brought outside the integral. The focal plane is divided into a region of interest

_{in}*Ω*and its complement

*Ω*. The region

^{c}*Ω*wholly contains and bounds the targeted beam shape

*I*.

_{t}#### 2.2. Normalization

It is helpful to cast the vector diffraction integrals into a form consisting of dimensionless variables by normalizing to the lateral extent of the input beam *I _{in}* and the targeted beam profile

*I*[21]. The Cartesian coordinates of the aperture plane (

_{t}*x*,

_{a}*y*) can be related to the

_{a}*x*- and

*y*-components of the wave vector by

where *R* is the radius of the entrance pupil, and *u* and *v* represent the normalized *k _{x}* and

*k*components of the wave vector, respectively. Likewise, the focal plane coordinates (

_{y}*x*,

_{f}*y*) are normalized by

_{f}*D*=

*mλ*(see Fig. 1), giving the transformed focal plane coordinates (

*ξ*,

*η*):

The size of *I _{t}* can be scaled conveniently by

*m*multiples of

*λ*. All free parameters of the system can be combined into the single variable,

*β*, given by [21]:

Substituting Eqs. (3) and (4) into Eq. (1) and algebraically manipulating gives the following normalized vectorial diffraction integrals for the field at the focal plane:

Where *q*(*u*, *ν*)=[1-(*NA*/*n*)^{2}(*u*
^{2}+*ν*
^{2})]^{1/2}, and *ω _{x}*=

*k*

_{max}

*Dξ*=

*βξ*and

*ω*=

_{y}*k*

_{max}

*Dη*=

*βη*have been introduced so the diffraction integrals in Eq. (5) can be evaluated as a Fourier transform.

#### 2.3. Chirp-z transform

In this work, the double integrals of Eq. (5) were evaluated using a two-dimensional (2D) chirp-z transform (CZT) with 512×512 sampling points in both the aperture and the focal planes, irrespective of the magnitude of *R* and *D*. Although 2D fast Fourier transform can be used, CZT is computationally faster and better suited for the present situation because it internalizes zero-padding and allows the spacing of sampling points in the aperture and focal planes to be set independently [22]. This greatly reduces the number of sampling points required when *R* and *D* differ substantially in magnitude, as in the present case.

Evaluating Eq. (5) effectively propagates the field *E _{x}*(

*u*,

*v*) forward into

*E*(

_{x}*ξ*,

*η*),

*E*(

_{y}*ξ*,

*η*), and

*E*(

_{z}*ξ*,

*η*). The fundamental operation of beam shaping involves applying constraints associated with

*I*to the field in the focal plane and then calculating a new DOE phase profile that comes closest to generating the reshaped field. A new DOE transmission profile is obtained by “backward propagating” to the pupil plane through an inverse-CZT applied to the reshaped field.

_{t}Because the forward and inverse CZTs are applied to a finite region of the DOE and focal planes, Gibbs artifacts are generated [18]. If these numerical errors are allowed to accumulate, they can degrade the uniformity of *I _{f}* or even cause the algorithm to diverge from a solution. Gibbs artifacts were suppressed by applying a Kaiser window to the amplitude of the focal field profiles immediately after they were computed using CZT [18].

#### 2.4. Method of generalized projections

For a given input beam *I _{in}* we seek a DOE phase function

*Φ*that generates a focal plane field distribution such that the sum of the

*x*-,

*y*-, and

*z*-polarized irradiance

*I*+

_{x}*I*+

_{y}*I*=

_{z}*I*, matches the targeted irradiance

_{f}*I*for all (

_{t}*ξ*,

*η*). An exact match is generally not possible because it is not known

*a priori*that an arbitrary

*I*is a solution to the wave equation [11]. This is the case for the present example because the targeted square irradiance profile requires a discontinuous drop in the field at the interface of

_{t}*Ω*and

*Ω*. The problem is further complicated because

^{c}*Φ*(

*u*,

*ν*) affects each of

*E*,

_{x}*E*, and

_{y}*E*, so the field components are not truly independent. Evaluating Eq. (5) and integrating

_{z}*I*,

_{x}*I*, and

_{y}*I*over the entire focal plane shows that their fractional power content is 0.74, 0.01 and 0.25, respectively, and these values are independent of

_{z}*Φ*(

*u*,

*ν*). So, high-

*NA*beam shaping demands that

*I*,

_{x}*I*, and

_{y}*I*are optimized collectively. This problem cannot be solved analytically, so iterative numerical techniques must be employed. The MGP is particularly well suited to the current problem because it can find solutions that closely satisfy sets of inconsistent and non-physical constraints [23]. In the MGP the optical field is repeatedly propagated forward and backward between the DOE and focal planes, and constraints associated with both domains are applied in each iteration until a satisfactory solution is found.

_{z}#### 2.5. Starting conditions

The goal is to design a phase-only DOE, so the initial transmission amplitude *T _{0}* is set to unity. The rate of convergence and quality of the solution can be greatly improved by initiating the vector diffraction algorithm using a well chosen starting DOE phase profile,

*Φ*

_{0}(

*u*,

*ν*). Geometrical optics based methods can be used to identify suitable starting DOEs [5]. Geometrical transformations have been applied successfully to obtain

*Φ*

_{0}(

*u*,

*ν*) analytically when

*I*and

_{in}*I*are either separable or axially symmetric [7]; however, the beam shaping example studied here is neither separable nor axially symmetric. To overcome this problem, a procedure described by Aagedal

_{t}*et al.*[7] was employed to obtain

*Φ*

_{0}(

*u*,

*ν*) as a combination of two separate DOEs that together achieve the required geometric beam transformation. The first element converts the axially symmetric

*I*into a standard Gaussian beam, which is separable and axially symmetric. The second element converts the Gaussian beam into a square-shaped super-Gaussian beam. The resulting

_{in}*Φ*

_{0}(

*u*,

*ν*) does not adequately reshape

*I*to

_{in}*I*under high-

_{t}*NA*, but it provides a good starting point for the vector diffraction algorithm.

#### 2.6. Algorithm flow

The analytically calculated starting DOE, *Φ _{0}*(

*u*,

*ν*), is substituted into Eq. (5). The diffraction integrals are then evaluated using the CZT giving |

*E*|exp(

_{x}*iϕ*), |

_{x}*E*|exp(

_{y}*iϕ*), and |

_{y}*E*|exp(

_{z}*iϕ*) in the focal plane, and the corresponding irradiance distribution is

_{z}*I*. The diffraction efficiency and uniformity error are calculated for the solution as

_{f}Under these definitions, a “perfect” solution would have *κ*=1 and *δ*=0.

Second, the constraint of the targeted beam shape *I _{t}* is applied. For the present example,

*I*is defined as

_{t}with the limits of *Ω* set to -0.25≤*ξ*,*η*≤0.25. Within the beam shaping area, total power is conserved and the irradiance is homogenized by setting the latter to its average value:

The *x*-polarized field within *Ω* is reshaped as

whereas the following are left unchanged: |*E _{x}*(

*ξ*,

*η*)| outside

*Ω*; |

*E*(

_{y}*ξ*,

*η*)| and |

*E*(

_{z}*ξ*,

*η*)| across all

*Ω*+

*Ω*; and the phases

^{c}*ϕ*(

*ξ*,

*η*) of all field components in

*Ω*+

*Ω*.

^{c}*γ*is an adjustable scalar that augments |

*E*(

_{x}*ξ*,

*η*)| in

*Ω*relative to that in

*Ω*. This operation provides a means for slowly pulling energy from

^{c}*Ω*into

^{c}*Ω*[8]. In this work

*γ*=1.03 was used for all iterations.

Third, an inverse CZT is applied to the reshaped *E _{x}*(

*ξ*,

*η*) to generate a complex DOE function

*T*

_{i+1}(

*u*,

*ν*)exp[

*iΦ*

_{i+1}(

*u*,

*ν*)]. Given that a phase-only DOE is required, we apply this constraint by resetting the transmission amplitude to

*T*while retaining the phase. The new complex DOE transmission function becomes

_{0}*T*(

_{0}*u*,

*ν*)exp[

*iΦ*

_{i+1}(

*u*,

*ν*)]. The electric field is then propagated forward again using the new DOE and the reshaped beam it generates is evaluated based on

*κ*and

*δ*. This process continues until the algorithm converges to a suitable solution or until a fixed number of iterations are completed.

## 3. Results and discussion

Equation (10) is intentionally configured so that the beam shaping constraint is only directly applied to the *x*-component of the field amplitude lying within *Ω*. This arrangement provides amplitude freedom outside the region of interest and phase freedom across the entire focal plane that help the algorithm converge to a solution [11, 24]. Additionally, it solves the problem of how to reshape three independent field components that are effectively coupled through a common DOE. The intended beam shape is applied repeatedly to the *x*-polarized field, as it contains the majority of the focused power, and only it is propagated backward to obtain the DOE phase function for the next iteration. The _{y}- and _{z}-polarized components of the focal field are reshaped indirectly when they are calculated by forward propagation through the new DOE. Repeated iterations effectively pull *E _{x}*(

*ξ*,

*η*),

*E*(

_{y}*ξ*,

*η*), and

*E*(

_{z}*ξ*,

*η*) toward distributions that collectively satisfy

*I*+

_{x}*I*+

_{y}*I*→

_{z}*I*.

_{t}We now discuss the results obtained when the vector diffraction algorithm was used to design a phase-only DOE that transforms a circularly apodized flat-top input beam into a focused square irradiance pattern of area 50*λ*×50*λ* (*D*=100*λ*). Figure 2 shows the normalized focal irradiance distributions of the three polarization components, *I _{x}*,

*I*, and

_{y}*I*, and the total focal irradiance distribution

_{z}*I*generated by the DOE phase profile of Fig. 3. This DOE and the associated irradiance distributions were obtained after 600 iterations. The overall beam shape is square as intended with

_{f}*δ*=7% and

*κ*=74.5%, indicating that it has good uniformity and power confinement within the region of interest.

In contrast, the irradiance distributions of the constituent polarizations are non-uniform. *I*
_{x} most resembles the targeted profile, but appears doubly concave, as though squeezed along the *x*-axis. Although *I _{x}* is non-zero across the coordinate axes,

*I*and

_{y}*I*have node(s) at these positions where their field amplitudes drop to zero.

_{z}*I*is most complex, appearing approximately four-fold symmetric with power concentrated in the corners of

_{y}*Ω*.

*I*exhibits two-fold symmetry with a single nodal plane lying along the

_{z}*y*-axis. The regions of high irradiance in

*I*and

_{y}*I*fill in around the edges of the

_{z}*x*-polarized profile making the total irradiance distribution

*I*uniform and square. These profiles show that the vector diffraction algorithm successfully generates a DOE for which all polarization components of the field are reshaped concurrently to achieve a targeted irradiance distribution under high-

_{f}*NA*focusing.

Figure 4 shows how the uniformity error and diffraction efficiency change during the calculation. The diffraction efficiency progressively increases because the parameter *γ*=1.03 causes power to transfer from *Ω ^{c}* into

*Ω*with each iteration. On the other hand, the uniformity error drops rapidly and reaches an apparent plateau after circa 500 iterations. It is known that high uniformity and high diffraction efficiency are mutually exclusive characteristics in beam shaping [19]. As a result, attempting to improve the diffraction efficiency beyond the level of 74% achieved at approximately 600 iterations caused the uniformity to erode. Obtaining solutions that are optimized in terms of both

*δ*and

*κ*could be achieved by extending the present vector diffraction algorithm through Tikhonov regularization theory [19]. It is noteworthy that the diffraction efficiency and uniformity are very poor for the first iteration. This results because the starting DOE was designed using a geometrical transformation method, which does not account for the vector character of the field. It underscores the importance then of using vector diffraction theory to achieve accurate beam shaping under high-

*NA*conditions, and it demonstrates the improvement that can be achieved in beam shaping using the present vectorial approach.

The problem of shaping a beam whose size approaches the diffraction limit was examined by repeating the calculations described above for *D*=50*λ*, 25*λ*, 10*λ*, and 5*λ*. The resulting focused irradiance distributions and the corresponding DOE phase profiles are shown in Fig. 5. Comparing the irradiance distributions reveals that the intended beam transformation can be achieved, even for targeted beam profiles having an edge-length of *D*/2=2.5*λ* (see Fig. 5(G)). The DOEs themselves have approximate four-fold symmetry with respect to rotation about the *z*-axis, as expected for a square target beam shape (consider also Fig. 3). The DOEs are also comprised of many concentric rings of steadily increasing phase reminiscent of a Fresnel lens. These concentric phase rings effectively negate some of the focusing power of the high-*NA* lens, so understandably their number and radial density decreases as the target beam size is reduced toward the diffraction limit.

As *D* decreases, the reshaped beam degrades in uniformity and sharpness at the boundary of *Ω*. The sharpness of the irradiance profiles was characterized empirically by fitting the normalized individual distributions to a super-Gaussian of the form

The order of the super-Gaussian, *N*, provides a measure of the sharpness of the beam profile at the boundary of the region of interest. Rather than applying Eq. (11) to all of *Ω*+*Ω ^{c}*, the fitting was restricted to a region within which

*I*exceeds 0.5. This procedure yields fits that more accurately describe the steepness of the profiles at the interface between

_{f}*Ω*and

*Ω*because it does not include irradiance fluctuations in

^{c}*Ω*that are necessarily part of any real solution to the wave equation. The beam shaping parameters

^{c}*κ*,

*δ*, and

*N*obtained for each value of

*D*are collated in Table 1. The data show that the beam uniformity, diffraction efficiency, and edge sharpness all deteriorate as the size of the reshaped beam is reduced toward the diffraction limit.

The results in Fig. 5 suggest that it becomes increasingly difficult to achieve a targeted irradiance distribution when the beam size becomes comparable to the diffraction limit, as has also been observed for scalar beam shaping [25]. This phenomenon can be understood as a manifestation of the uncertainty principle [21]. The irradiance distribution *I _{t}* that can be achieved is inherently limited by the finite spatial bandwidth of the input beam

*I*and the limited range of wave vectors over which focusing occurs, as quantified by the

_{in}*NA*. The limit in the achievable beam shape can be expressed as [21]

where

The particular usefulness of *β* now becomes apparent. Given that *β* is determined by all free parameters of the system (*R*, *D*, *λ*, and *NA*), it provides a single measure of the difficulty of the beam shaping problem. Values of *β ^{2}ΔI_{in}ΔI_{t}* for the beam profiles of Fig. 5 are also included in Table 1. Because

*β*depends quadratically on

^{2}ΔI_{in}ΔI_{t}*D*, it drops rapidly within this series and is most comparable to unity at

*D*=5

*λ*, for which the reshaped beam quality is poorest. These data and Eq. (12) imply then that

*I*will differ increasingly from

_{f}*I*as the targeted beam size is decreased toward the diffraction limited spot size, with all other parameters kept fixed.

_{t}## 4. Conclusion

A vector diffraction algorithm was developed for designing phase-only DOEs that reshape beams focused under high-*NA* conditions. The algorithm accounts for depolarization effects that occur under high-*NA* focusing by relating the DOE complex transmittance function and the electric field in the focal plane using the vector diffraction integrals. The algorithm was applied in the design of a phase-only DOE that reshapes the focused irradiance distribution of a circularly apodized flat-top input beam into a uniform square profile when focused using a 1.4-*NA* objective lens. We observe that beam uniformity, diffraction efficiency, and edge sharpness all degrade as the size of the targeted flat-top beam is decreased. This suggests that beam shaping becomes increasingly difficult as the area of the targeted irradiance distribution approaches that of the diffraction limited spot size. There are many possibilities for extending the method reported here. A wider range of focused beam shapes could be considered by appropriately modifying the constraints used to define the targeted irradiance distribution. The search for solutions could be made more general by including other free parameters. For example, one could allow for freedom in the size of the homogenized area, so that the targeted size is optimized along with diffraction efficiency and uniformity. Three-dimensional beam shaping could be achieved by applying constraints in multiple planes. This work may also be useful for extending methods employed in the design of phase masks for high-resolution photo-lithography.

## Acknowledgements

We wish to thank Marcel Leutenegger for helpful discussions pertaining to this work. We also thank the reviewers for their many helpful and insightful comments.

## References and links

**1. **C. Dorrer and J. D. Zuegel, “Design and analysis of binary beam shapers using error diffusion,” J. Opt. Soc. Am. B **24**, 1268–1275 (2007). [CrossRef]

**2. **V. Ramírez-Sánchez and G. Piquero, “Global beam shaping with nonuniformily polarized beams: a proposal,” Appl. Opt. **45**, 8902–8906 (2006). [CrossRef] [PubMed]

**3. **L. A. Romero and F. M. Dickey, “Lossless laser beam shaping,” J. Opt. Soc. Am. A **13**, 751–760 (1995). [CrossRef]

**4. **F. M. Dickey, S. C. Holswade, and D. L. Shealy, *Laser Beam Shaping Applications* (Taylor & Francis Group, Boca Raton, 2006).

**5. **O. Bryngdhal, “Geometrical transformations in optics,” J. Opt. Soc. Am. **64**, 1092–1099 (1974). [CrossRef]

**6. **D. L. Shealy and S. H. Chao, “Geometric optics-based design of laser beam shapers,” Opt. Eng. **42**, 3123–3138 (2003). [CrossRef]

**7. **H. Aagedal, M. Schmid, S. Egner, J. Müller-Quade, T. Beth, and F. Wyrowski, “Analytical beam shaping with application to laser-diode arrays,” J. Opt. Soc. Am. A **14**, 1549–1553 (1997). [CrossRef]

**8. **M. Johansson and J. Bengtsson, “Robust design method for highly efficient beam-shaping diffractive optical elements using iterative-Fourier-transform algorithm with soft operations,” J. Mod. Opt. **47**, 1385–1398 (2000). [CrossRef]

**9. **V. V. Kotlyar, P. G. Seraphimovich, and V. A. Soifer, “An iterative algorithm for designing diffractive optical elements with regularization,” Opt. Laser Eng. **29**, 261–268 (1998). [CrossRef]

**10. **J. S. Liu and M. R. Taghizadeh, “Iterative algorithm for the design of diffractive phase elements for laser beam shaping,” Opt. Lett. **27**, 1463–1465 (2002). [CrossRef]

**11. **F. Wyrowski, “Diffractive optical elements: iterative calculations of quantized, blazed phase structures,” J. Opt. Soc. Am. A **7**, 961–969 (1990). [CrossRef]

**12. **G. Zhou, X. Yuan, P. Dowd, Y. L. Lam, and Y. C. Chan, “Design of diffractive phase elements for beam shaping: hybrid approach,” J. Opt. Soc. Am. A **18**, 791–800 (2001). [CrossRef]

**13. **M. Kuittinen, P. Vahimaa, M. Honkanen, and J. Turunen, “Beam shaping in the nonparaxial domain of diffractive optics,” Appl. Opt. **36**, 2034–2041 (1997). [CrossRef] [PubMed]

**14. **T. G. Jabbour and S. M. Kuebler, “Vector diffraction analysis of high numerical aperture focused beams modified by two- and three-zone annular multi-phase plates,” Opt. Express **14**, 1033–1043 (2006). [CrossRef] [PubMed]

**15. **Y. Zhao, Y.-P. Li, and Q.-G. Zhou, “Vector iterative algorithm for the design of diffractive optical elements applied to uniform illumination,” Opt. Lett. **29**, 664–666 (2004). [CrossRef] [PubMed]

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

**17. **A. Levi and H. Stark, “Image restoration by the method of generalized projections with application to restoration from magnitude,” J. Opt. Soc. Am. A **1**, 932–943 (1984). [CrossRef]

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

**19. **H. Kim, B. Yang, and B. Lee, “Iterative Fourier transform algorithm with regularization for the optimal design of diffractive optical elements,” J. Opt. Soc. Am. A **21**, 2353–2365 (2004). [CrossRef]

**20. **R. Kant, “An analytical solution of vector diffraction for focusing optical systems with Seidal aberrations I. Spherical aberration, curvature of field, and distortion,” J. Mod. Opt. **40**, 2293–2310 (1993). [CrossRef]

**21. **F. M. Dickey and D. L. Shealy, *Laser Beam Shaping* (Marcel Dekker, Inc., New York, 2000). [CrossRef]

**22. **M. Leutenegger, R. Rao, R. A. Leitgeh, and T. Lasser, “Fast focus field calculations,” Opt. Express **14**, 11277–11291 (2006). [CrossRef] [PubMed]

**23. **R. Piestun and J. Shamir, “Synthesis of three-dimensional light fields and applications,” Proceedings of the IEEE **90**, 222–244 (2002). [CrossRef]

**24. **F. Wyrowski, “Diffraction efficiency of analog and quantized digital amplitude holograms: analysis and manipulation,” J. Opt. Soc. Am. A **7**, 383–393 (1990). [CrossRef]

**25. **J. S. Liu, A. J. Caley, and M. R. Taghizadeh, “Symmetrical iterative Fourier-transform algorithm using both phase and amplitude freedoms,” Opt. Commun. **267**, 347–355 (2006). [CrossRef]