Open Access
Add to BibSonomyAbstract
We have designed high-efficiency finite-aperture diffractive optical elements (DOE’s) with features on the order of or smaller than the wavelength of the incident illumination. The use of scalar diffraction theory is generally not considered valid for the design of DOE’s with such features. However, we have found several cases in which the use of a scalar-based design is, in fact, quite accurate. We also present a modified scalar-based iterative design method that incorporates the angular spectrum approach to design diffractive optical elements that operate in the near-field and have sub-wavelength features. We call this design method the iterative angular spectrum approach (IASA). Upon comparison with a rigorous electromagnetic analysis technique, specifically, the finite difference time-domain method (FDTD), we find that our scalar-based design method is surprisingly valid for DOE’s having sub-wavelength features.
©2001 Optical Society of America
1. Introduction
Scalar diffraction theory is often used in the design and analysis of diffractive optical elements (DOE’s). It is straightforward to use and the computational burden is minimal. Generally, the use of scalar diffraction theory is only considered valid if the DOE minimum feature size is much greater than the wavelength of incident light [1]. However, we have found several applications in which a scalar-based diffraction theory is surprisingly accurate, even though the minimum features are smaller than the wavelength of the incident light.
In this paper, we discuss design examples of beamfanners that illustrate the use of a scalar-based design algorithm and explore limitations of the use of scalar-based diffraction theory to design DOE’s. There are two fundamental assumptions associated with the application of scalar diffraction theory to design DOE’s. The first is that the optical field just past the DOE can be related to the incident field by a simple transmission function [2]. The second involves the choice of method to propagate the field to the plane of interest. In most applications presented in the literature, the plane of interest is located in the far-field. In our case the plane of interest is in the near field so the usual selection of Fraunhofer diffraction for the propagation method is not accurate. Fresnel propagation is commonly used in the near-field. This method, however, is still only approximate. We therefore use the angular spectrum approach [2,3], which is rigorous within scalar diffraction theory and is formally equivalent to electromagnetic TE electric field and TM magnetic field propagation. To make the distinction clear, we use the term angular spectrum (AS) scalar diffraction theory to emphasize that the angular spectrum approach is used as the method of propagation.
The limits of scalar diffraction theory as the DOE minimum feature size approaches the wavelength of illuminating light have not been fully quantified in the literature, and only a few authors have attempted to assess scalar theory in this regime. For example, Gremaux and Gallagher investigated the limits of scalar diffraction theory for perfectly conducting gratings [4] while Pommet et. al [5], examined the limits of scalar diffraction theory as applied to phase gratings. In the case of Pommet et. al, they compared the diffraction characteristics of single-level and multilevel dielectric gratings predicted by scalar transmission theory with those obtained by rigorous coupled-wave analysis (RCWA) [6], and found that, in general, the error of scalar theory is significant when the feature size is less than 14 wavelengths. They also report that the error is minimized when the fill factor approaches 50% for DOE’s with smaller features of 2 wavelengths. In such cases, the larger period of the grating replaces the feature size as the condition of validity for scalar diffraction theory.
In this paper, we consider the case of phase-only finite aperture DOE’s with subwavelength minimum features. To aid in this analysis, we introduce the concept of zones in which we define a zone as a region in which the DOE profile is quasi-continuous, i.e., the phase levels of adjacent subwavelength minimum features within a zone do not change abruptly. (Note that this is different from the concept of Fresnel zones, which are related to the 2π phase resets in a Fresnel lens.) A key finding of this paper is that AS scalar theory is reasonably accurate even for subwavelength minimum features provided that the zone size of the DOE is at least on the order of the free-space wavelength.
Based on this result, we present a bidirectional modified iterative design algorithm similar, in principle, to the Gerchberg-Saxton algorithm [7], the iterative Fresnel-transform algorithm (IFTA), and the iterative Fourier-transform algorithm [8–12]. Adopting the nomenclature used by Mait [13], bidirectional means that updating the DOE profile requires knowledge of the inversion of the operation used to obtain the field the DOE produces and how variations in the response affect the DOE. In contrast, unidirectional algorithms such as simulated annealing [14–18] and micro-genetic algorithms [19], characterize the DOE by a finite set of quantized parameters such as phase levels, and do not require an inversion operation to update the DOE.
Our design method incorporates the angular spectrum approach in a bidirectional, iterative algorithm. We call this method the iterative angular spectrum approach (IASA). Since the angular spectrum approach requires no approximations (as do Fresnel and Fraunhofer propagation), the main approximation in our design algorithm is the use of a transmission function to obtain the optical field immediately past the DOE. The use of a transmission function permits the back-propagated field to be related to the DOE phase profile, which enables the algorithm to be bidirectional. The success or failure of our design approach for subwavelength minimum features is therefore primarily related to the validity of the transmission function approximation, and not the propagation method. The results presented in Sect. 3.3 for 1-to-3 and 1-to-4 beamfanners indicate that this approximation can be surprisingly accurate, even for subwavelength minimum features.
In Sec. 2, we present designs of 1–2 beamfanners and assess the accuracy of a scalar diffraction approach used to heuristically design and analyze the DOE’s by comparison with rigorous results. In Sec. 3, we present our AS scalar-based iterative design method, IASA, design examples, and a comparison with results from a rigorous electromagnetic analysis. In Sec. 4, we provide a summary and concluding remarks.
2. Heuristic design and analysis of 1–2 beamfanners
Consider a DOE that functions as a 1–2 beamfanner as shown in Fig. 1. The DOE’s function is to maximize the amount of light directed into the two apertures in what we will refer to as the observation plane. Note that the observation plane is in the near-field. A plane-wave with a free-space wavelength, λ0, of 5 µm is normally incident on a silicon-air interface bounded by a 50 µm (10 λ0) metallic aperture having infinite conductivity. The minimum feature size of the DOE is chosen to be 1 µm (0.2 λ0) and the observation plane is located at a distance of 100 µm in air (20 λ0). Next, we consider a heuristic method to design a DOE to perform the 1–2 beam-fanning function.
To begin the design, we take an infinite grating with spatial period, Λ, and truncate it to fit the bounds of the finite aperture. The fill factor of the grating is set to 50% and its depth, d, is set such that the power of the ±1 diffracted orders are maximized according to scalar theory. This occurs when the relative phase shift between adjacent zones is π radians, which corresponds to d=λ0/(2Δn) for normally incident light and Δn is the difference in refractive index between air and silicon. Here, we define a zone as a region in which the DOE profile is locally continuous, e.g., for a binary grating with a 50% fill factor, a zone would constitute half the grating period. Next, a quadratic phase profile is added to the grating to focus the light in the near field observation plane.
For scalar-based analysis of DOE diffraction, we assume a transmission function of the form t(x)=τ exp(jφ(x)) rect(x/L) in which L is the width of the finite aperture and rect(x/L) equals unity for |x|<L/2 and is zero otherwise. The profile phase, φ(x), can be expressed as φ(x)=2πΔnd(x)/λ0 in which d(x) is the etch depth along the DOE. The Fresnel transmission coefficient, τ, is given by τ=2n1cos(θ1)/(n1 cos(θ1)+n2 cos(θ2)) in which n1 and n2 are the refractive indices of silicon and air, respectively, θ1 is the angle of the incidence with respect to the optical axis and θ2 is the refracted angle as calculated by Snell’s law. In every case in this paper, we only consider normally incident light, for which the Fresnel transmission coefficient reduces to τ=2 n1/(n1+n2). The field incident on a DOE is multiplied by the transmission function, t(x), which gives the field just past the DOE interface. This field is then propagated to a chosen observation plane using the angular spectrum approach. Once the field is known in this observation plane, all other field component can be obtained via Maxwell’s equations, and the power and the diffraction efficiency can be calculated. By using AS scalar-based theory, the goal is to determine what DOE profile, d(x), best performs its intended function. In the heuristic designs presented, we specifically design 1–2 beamfanners.
Given that the DOE’s have sub-wavelength features, it may appear counterintuitive to use a scalar-based approach in their design. It turns out, however, that several of these designs are very accurate under rigorous scrutiny. All scalar-based results were compared to those obtained using a rigorous finite-difference time-domain method (FDTD) for both TE and TM incident illumination [20]. The FDTD analysis, using Mur 2nd order boundary conditions and the total-field/scattered-field formulation to introduce the incident fields, employed a time-marching scheme to analyze the electric and magnetic field values at sampled points within the FDTD grid, which enclosed the entire DOE structure. At a plane immediately past the DOE in the exiting medium intersecting the FDTD grid, the field was propagated to the observation plane using the electromagnetic angular spectrum approach [2,3].
In the test cases, we examined the effect of varying the spatial period of the grating (and hence the zone size) on the accuracy of applying AS scalar diffraction to analyze DOE’s. The values of Λ chosen were 8 λ0, 4 λ0, 2 λ0, 1.6 λ0, 1.4 λ0, and 1.2 λ0, and the resultant DOE profiles for these cases are shown in Fig. 2.
Fig. 3 shows the corresponding near-field irradiance in the observation plane. Given that the distance of this observation plane is 100 µm from the DOE interface and the location of the diffracted orders in that plane, note that the periods are self-consistent with those predicted by the grating equation: λ=Λsin(θ) [1], in which λ is the reduced wavelength in silicon (λ0/n), Λ is the grating period, and θ is the angle of the first diffracted order with respect to the optical axis, i.e., θ=tan-1(s/2z), s is the separation between the diffracted orders in the observation plane, and z is the distance to that plane. In Fig. 3, note that for larger grating periods the AS scalar results agree very well with those predicted by FDTD, whereas for smaller grating periods the validity of AS scalar theory appears to break down.
To quantify these results, we define the diffraction efficiency, η, to be the fraction of the incident light that is diffracted into the photodetector apertures in the observation plane, that is
in which I is the irradiance given by I=|Re{}·n̂| and is the complex Poynting vector in the observation plane and is the outward unit normal from the observation plane. The width of each of the photodetector apertures is chosen to be 15 µm. We also define the error in terms of diffraction efficiency between rigorous TE or TM results with respect to AS scalar results as
and the difference between the TM result with respect to the TE result as
Table 1 lists the diffraction efficiencies calculated from AS scalar theory and from FDTD TE and TM cases as well as the errors.
Table 1. Diffraction efficiencies and errors for heuristically designed 1–2 beamfanners
For smaller grating periods, the deviation of AS scalar results from those predicted by FDTD becomes more significant. This may indicate that AS scalar theory works well provided that the spatial period, Λ, is not much less than twice the free-space wavelength, i.e., Λ≥2λ0. This assessment appears to indicate that perhaps the validity and accuracy of AS scalar diffraction theory is slightly broader than other results reported in the literature [5]. In terms of zones, AS scalar diffraction theory is accurate, i.e., less than 10% error, for zone sizes greater than or equal to the free space wavelength of the incident light. Also note that the sub-wavelength minimum feature size is the same for each beamfanner in the analysis. Therefore, one may conclude that the validity of AS scalar diffraction theory does not necessarily depend on DOE minimum feature size but rather the size of the DOE zones.
Like AS scalar diffraction, our implementation of the FDTD method also incorporates the rigorous angular spectrum approach to propagate the field just past the DOE to the observation plane. Therefore, we conclude that it is the assumption of a simple transmission function that causes the breakdown of the validity of AS scalar theory in comparison with FDTD. We therefore examine the fields immediately past the DOE structure in the exit medium and compare the field amplitudes and phases predicted by AS scalar theory with those predicted by FDTD. Fig. 4 and Fig. 5 show the results of the comparison for the phases and amplitudes, respectively, for all the test cases considered. Note that the phases calculated by FDTD agree very well with AS scalar theory for large grating periods and begin to deviate to a certain extent as the grating period progressively decreases relative to λ0. The field amplitudes calculated via FDTD, however, differ more dramatically from the AS scalar result, which is simply a rectangle function scaled by the Fresnel transmission coefficient. This apparent disagreement in field magnitudes can be reconciled somewhat by examining only those field components that actually propagate to the observation plane, i.e., which are not evanescent, as will be shown next.
First consider the angular spectra of the fields at a plane just past the DOE, which are shown in Fig. 6. Note that for the cases of Λ equaling 8 λ0, 4 λ0, and 2 λ0, the magnitudes of the angular spectra calculated from rigorous analysis agree well with what was calculated from the scalar analysis while the cases 1.6 λ0, 1.4 λ0, and 1.2 λ0, tend to be somewhat different. Also note that the spatial frequency components outside the range of [-1/λ0,1/λ0] are evanescent. These limits are denoted by vertical red lines.
From the angular spectrum calculations, we can examine only those field components just past the DOE that actually propagate to the observation plane. Only the propagating field contributes to the irradiance in the observation plane. As shown in Fig. 7, examination of the propagating field amplitudes just past the DOE reveals that scalar and rigorous results are in much better agreement than the total fields of Fig. 4.
As shown in Fig. 8, comparison of the phases for the propagating field just past the DOE is not much different from the total field phase case. The analysis of the propagating field amplitudes just past the DOE indicate that the scalar and rigorous results are in reasonable agreement for the cases of Λ equaling 8 λ0, 4 λ0, and 2 λ0. However, AS scalar results tend to differ significantly from rigorous results when the grating period is less than 2 λ0, which is consistent with what we observe in the near-field diffraction patterns. From this, we conclude that the amplitude and the phase of the propagating field just past the DOE is the most dominant factor in determining whether AS scalar diffraction theory is valid in our test cases. This also suggests that any DOE profile that generates the desired propagating field just past the DOE, regardless of the evanescent field content of the total field, will yield the same diffracted beam (with the possible exception that the overall diffraction efficiency may be different). Hence the evanescent field components represent a general design freedom for creating a diffractive optical element.
3. Iterative angular spectrum algorithm
Having examined the validity of AS scalar diffraction theory, we present a scalar-based design algorithm for DOE’s with sub-wavelength minimum feature sizes.
3.1 Design method
An iterative angular spectrum approach (IASA) which is a modification of the Gerchberg-Saxton algorithm [7], or the Iterative Fourier Transform Algorithm (IFTA) [8–12], is used in the design of the beamfanners discussed below. The steps to IASA are as follows:
1. Assume initial DOE etch depth profile.
2. Calculate transmission function and then the field just past DOE.
3. Obtain angular spectrum of field via Fast-Fourier Transform (FFT).
4. Forward propagate angular spectrum.
5. Determine field in the observation plane via inverse FFT.
6. Compare to desired diffraction pattern.
7. If diffraction pattern is satisfactory, IASA is complete, otherwise modify field amplitude according to predetermined weight functions.
8. Obtain angular spectrum of this field.
9. Back-propagate angular spectrum to DOE object plane.
10. Add evanescent components from Step 3 to obtain total angular spectrum.
11. Fourier transform spectrum to obtain field just past DOE.
12. Modify DOE profile applying appropriate constraints (i.e. phase-only DOE bounded by a finite aperture).
13. Go back to step 2.
An interesting aspect of IASA is contained in Step 10 in which the evanescent angular spectrum components of the previous iteration are added to the back-propagated angular spectrum (which of course has no evanescent components). Since subwavelength features strongly affect the evanescent angular spectrum components, our motivation for this step is to recover subwavelength feature information. An alternative approach is to sample the phase of the field in Step 12 such that the phase is constant over each subwavelength feature. The outcome is a stairstep approximation to the phase of the field, with the step size being the same as the minimum feature size. Initial results with this approach indicate similar performance to the algorithm as formulated above. However, we have not fully explored this issue and view it as a fruitful area for future research.
3.2 1–2 beamfanner design
In the 1–2 beamfanner design, we assume that a unit amplitude plane wave is normally incident on the finite aperture DOE. The DOE splits the light between two photodetectors in an observation plane in the near field. The wavelength of light considered is 5 µm, in free space, and the silicon substrate has minimum features of 1 µm bounded by a 50 µm aperture. We assume that the etch depth levels are not quantized, although they could be in a post-processing step. The apertures of each photodetector in the observation plane are chosen to have a width of 15 µm. To assess the range in which the use of IASA is valid and accurate, we designed beamfanners that would require grating periods comparable in size and smaller than the wavelength in free space. The cases presented below are beamfanners intended to split light with peak separations of 25, 50, 100, 150, 200, and 250 µm in an observation plane located 340 µm in the silicon from the air-silicon interface. As will be seen below, the choice of these geometrical parameters yields beamfanner designs similar to the heuristic designs of Section 2. The motivation to propagate from air to silicon (rather than silicon to air as before) is to illustrate that the same general conclusions are valid, regardless of which side of the DOE interface constitutes the high index material..
IASA was used to design finite aperture DOE’s that split light in the plane of interest with the peak separations mentioned above. The assumed initial profile was a flat interface bounded by a 50 µm finite aperture in an infinitely conducting plane. [Other initial profiles (such as random and quadratic phase profiles) yielded similar results to those given below.] Magnitude constraints are imposed through the use of weighting functions so that the field amplitudes are increased each iteration at positions in the observation plane corresponding to the locations of the photodetector apertures. Specifically, the field in the detector apertures is multiplied by (1+ε) while the field outside the apertures is weighted by (1-ε), where ε≪1. We also found that we had to impose an additional constraint to prevent the algorithm from converging to two adjacent microlenses that focus light independently of one another. To prevent this outcome (which was important for our eventual application), only one half (either left or right) of the DOE is illuminated and the energy is equalized about the central position (x=0) in the observation plane. All of the constraints are applied in each iteration until an acceptable irradiance pattern is obtained. For all cases considered, IASA produced an acceptable irradiance pattern within 1000 iterations.
The resultant DOE etch depth profiles produced by IASA that yield the desired peak separations in the observation plane are shown in Fig. 9. Note that each etch depth profile looks progressively more like a grating (with possibly several 2π phase resets) with a quadratic envelope in which the grating performs the beamfanning function while the quadratic phase component focuses light in the observation plane. The etch depths in each case between each zone are roughly a π phase difference, corresponding to an etch depth between adjacent peaks and troughs of 1.03 µm in these cases, which is an expected result for a binary grating having the maximum diffraction efficiency [5]. Given the location of the observation plane and the peaks of the diffracted orders, the periods are consistent with those predicted by the grating equation. Upon inspection of the DOE profiles, the ratios of the grating periods to the free-space wavelength are approximately 8, 4, 2, 1.4, 1.04, and 0.85, corresponding to peak separations of 25, 50, 100, 150, 200, and 250 µm, respectively.
Again, a rigorous electromagnetic analysis based on the FDTD method was performed for both TE and TM cases to assess the accuracy of AS scalar diffraction theory. Fig. 10 shows a comparison of the irradiances for the AS scalar-predicted result and TE and TM results obtained via FDTD. Note that for the closer peak separations, 25, 50, and 100 µm, the AS scalar results from IASA agree very well with those predicted by rigorous electromagnetic theory. For the larger peak separations of 150, 200, and 250 µm, deviation of IASA results from those predicted by FDTD becomes progressively more significant. The diffraction efficiencies and errors, as defined previously, are shown in Table 2. Note that the errors of TE and TM with respect to scalar are less than 6.5 percent for cases where the peak separations are 25, 50, and 100µm. In these cases, IASA accurately designs the DOE’s for these functions. However, for larger peak separations (and hence smaller zone sizes) the errors become significant, indicating that the validity of AS scalar theory breaks down in these cases. Also, the diffraction efficiency predicted by AS scalar theory is much higher than what is calculated by FDTD for these larger peak separations.
Table 2. Diffraction efficiencies and errors for 1-2 beamfanners designed via IASA
3.3 1–3 and 1–4 beamfanner designs
The advantage of using IASA over a heuristic method is that the former does not require prior knowledge of the form of the DOE profile. Although 1–2 beamfanner designs can be done by either heuristic means or numerically with IASA, this application falls into a very limited class of diffractive optic design problems. For example, designs for a 1–3 and a 1–4 beamfanner cannot be done heuristically, and must be done numerically. Note that IASA, in general, could also be used for higher-order 1-N beamfanner designs in which N is the number of diffracted peaks, but in this paper we restrict ourselves to the following 1–3 and 1–4 beamfanner designs.
Consider an air-silicon interface, as before, in which silicon is the exiting medium and consider an observation plane located 340 µm from the interface in the silicon. Let the finite aperture remain 50 µm wide and consider 1 µm minimum feature sizes. For the designs of 1–3 and 1–4 beamfanners, the weighting functions are very similar to the 1–2 beamfanner designs with the difference that the field amplitudes are equally maximized at more locations. Fig. 11 shows the resultant DOE profiles with IASA and the corresponding near-field diffraction patterns in the observation plane. Table 3 lists the diffraction efficiencies calculated from AS scalar theory and for FDTD TE and TM cases as well as the previously defined errors. It appears that the use of IASA is, in fact, accurate since there is excellent agreement between the scalar and rigorous results.
Table 3. Diffraction efficiencies and errors for 1–3 and 1–4 beamfanners
4. Summary
In summary, it appears that AS scalar diffraction theory is applicable in several cases for finite aperture DOE design with features on the order of or less than the wavelength of illuminating light. Regardless of whether the exiting medium has a high or low refractive index, it appears that the AS scalar model is a viable method provided that the zone size is not less than the free space wavelength. We also presented an iterative design method, namely IASA, which is surprisingly accurate for designing finite aperture DOE’s with sub-wavelength features. Also, IASA can be used to design finite aperture DOE’s which cannot be designed heuristically, such as 1–3 and 1–4 or higher-order beamfanners.
We examined the phases and amplitudes of the total fields and propagating fields immediately after the finite aperture DOE structures. Future investigations would include more extensive analysis of the fields near DOE interfaces to more fully assess the limits of scalar-based diffraction theory. This investigation would also include examining the effects of DOE zone edges on the fields just past a DOE and may further explain the reasons why as well as the extent to which AS scalar diffraction theory is valid.
Acknowledgments
This work was supported by National Science Foundation CAREER Award ECS-9625040 and grant EPS-9720653.
References and links
1. M. Born and E. Wolf, Principles of Optics, 3rd ed. (Pergamon Press, New York, 1965).
2. J. Goodman, Introduction to Fourier Optics (McGraw-Hill., New York, 1968).
3. G.S. Smith. An Introduction to Classical Electromagnetic Radiation (Cambridge University Press, Cambridge, 1997).
4. D. A. Gremaux and N. C. Gallagher, “Limits of scalar diffraction theory for conducting gratings,” Appl. Opt. 32, 1948–1953 (1993). [CrossRef] [PubMed]
5. D.A. Pommet, M.G. Moharam, and E.B. Grann, “Limits of scalar diffraction theory for diffractive phase elements”, Opt. Lett. 11, 1827–1834 (1995).
6. M. G. Moharam and T. K. Gaylord, “Diffraction analysis of surface-relief gratings,” J. Opt. Soc. Am. A 72, 1385–1392 (1982). [CrossRef]
7. R.W. Gerchberg and W.O. Saxton. “A practical algorithm for the determination of phase from image and diffraction plane pictures” Optik 35, 237–246 (1971).
8. Pierre St. Hilaire, “Phase profiles for holographic stereograms,” Opt. Eng. 34, 83–89 (1995). [CrossRef]
9. N.C. Gallagher and B. Liu, “Method for computing kinoforms that reduces image reconstruction error,” Appl. Opt. 12, 2328–2335 (1973). [CrossRef] [PubMed]
10. J.R. Fienup, “Iterative method applied to image reconstruction and to computer-generated holograms,” Opt. Eng. 19, 297–306 (1980).
11. F. Wyrowski, “Iterative Fourier-transform algorithm applied to computer holography,” J. Opt. Soc. Am. A 5, 1058–1065 (1988). [CrossRef]
12. F. Wyrowski and O. Bryngdahl, “Digital holography as part of diffractive optics,” Rep. Prog. Phys. 54, 1481–1571 (1991). [CrossRef]
13. J. N. Mait, “Understanding diffractive optic design in the scalar domain,” J. Opt. Soc. Am. A 12, 2145–2158 (1995). [CrossRef]
14. I.O. Bohachevsk;y, M.E. Johnson, and M.L. Stein, “Generalized simulated annealing for function optimization,” Technometrics 28, 209–217 (1986). [CrossRef]
15. Y. Lin, T.J. Kessler, and G.N. Lawrence, “Design of continuous surface-relief phase plates by surface-based simulated annealing to achieve control of focal-plane irradiance,” Opt. Lett. 21, 1703–1705 (1996). [CrossRef] [PubMed]
16. B. K. Jennison, J. P. Allebach, and D. W. Sweeney, “Iterative approaches to computer-generated holography,” Opt. Eng. 28, 629–637 (1989).
17. J. Turunen, A. Vasara, and J. Westerholm, “Kinoform phase relief synthesis: a stochastic method,” Opt. Eng. 28, 1162–1167 (1989).
18. M.R. Feldman and C.C. Guest, “High-efficiency hologram encoding for generation of spot arrays,” Opt. Lett. 14, 479–481 (1989). [CrossRef] [PubMed]
19. J. Jiang and G. Nordin, “A rigorous unidirectional method for designing finite aperture diffractive optical elements,” Optical Express , 7, 237–242 (2000), http://www.opticsexpress.org/oearchive/source/23164.htm. [CrossRef]
20. A. Taflove, Computational Electrodynamics: The Finite-Difference Time-Domain Method (Artech House, Boston, 1995).
