Numerical Evaluation of Diffraction Integrals.

This paper describes a simple numerical integration method for diffraction integrals which is based on elementary geometrical considerations of the manner in which different portions of the incident wavefront contribute to the diffracted field. The method is applicable in a wide range of cases as the assumptions regarding the type of integral are minimal, and the results are accurate even when the wavefront is divided into only a relatively small number of summation elements. Higher accuracies can be achieved by increasing the number of summation elements and/or incorporating Simpson's rule into the basic integration formula. The use of the method is illustrated by numerical examples based on Fresnel's diffraction integrals for circular apertures and apertures bounded by infinite straight lines (slits, half planes). In the latter cases, the numerical integration formula is reduced to a simple recursion formula, so that there is no need to perform repetitive summations for every point of the diffraction profile.

1. Basic Equations

Diffraction problems that can be described in two dimensions as indicated in Fig. 1 usually lead to complex integrals of the type where A and B account for the amplitude and phase of the optical field at a point of observation P = (x,z) due to a point Q = (ξ,0) located inside a diffracting aperture. In practice, the phase term B oscillates rapidly inside the range of integration while the amplitude term A varies slowly. The phase term B is often, but not always, a sinusoidal function of the form ei(), where k = 2π/λ is the circular wavenumber of monochromatic light with wavelength λ.

Fig. 1

Cross section of a plane diffracting aperture (xz-plane).

When closed analytical expressions for U are not available it is sometimes possible to find approximate solutions by the method of stationary phase [1,2]. However, in this author’s experience, the results obtained can be unreliable. A potential source of error is the basic premise of the method itself; namely, that the rapid oscillations of the phase term B nullify each other except in the vicinity of stationary points. This is presumed true on account of the slow variation of the amplitude term A, but inconsistent with the well-known fact that the corresponding series encountered in connection with Fresnel’s zone construction, A1−A2 + A3 −…− A, has a finite value even when the terms alter their absolute values very gradually [3,4]. Accordingly, substantial errors can occur if the contribution of the stationary points is weak. As shown in Appendix A, the stationary-phase method fails completely when applied to Fresnel diffraction at a slit or half plane.

For these reasons, it may be preferable to use numerical integration. In doing so, the width of the summation elements, ∆ξ = h, must be sufficiently small to ensure that all oscillations of the phase term B are accurately sampled. On the other hand, the computations would be unnecessarily complicated if h is too small. A general rule for choosing h can be established as follows. As we are dealing with diffraction, the period of the oscillations does not depend on the specific functional form of B but only on the path length P0Q + QP, where P0 = (x0,z0) is the source point shown in Fig. 1. When the point Q is moved along the x-axis by an increment h, this path length changes by where it is assumed that (ξ−x0)2 ≪ z02, (ξ−x)2≪z2, and terms in h2 are ignored. Hence it follows from the quarter-wave criterion that reliable results can be expected when h is chosen so that δ<λ/4.

We now divide the aperture halfwidth a into N summation elements bounded by equidistant points Q, as illustrated in Figs. 2a, b, and 4, and define where O is the coordinate origin, Q0 is the projection of the point of observation P onto the x-axis, and m, n, and N are integers. Accordingly, Eq. (1) can be replaced by the quadrature formula where are the values of A and B at the midpoints of the summation elements and the limits of summation depend on the diffraction problem being considered. The value of N to be used in these formulae can be estimated from Eq. (2) by assuming a distant source (|z0|≫z) so that δ ≈ h(ξ−x)/z = a2n/(N2z). If this is to be less than λ/4 for every summation element used for the computations, a good upper limit for n is 3N,1 and then one finds where u = ka2/z is the familiar configuration parameter of Fresnel’s diffraction theory for |z0|≫z. As it is well known that diffraction patterns pertaining to large values of u are highly structured [4], this result makes good sense in that it stipulates narrower summation elements when u increases. The accuracy of the numerical computations can also be improved by replacing the value of B in Eq. (3c) with the value corresponding to Simpson’s rule, By means of trial computations, it was found that this substitution can result in a tenfold improvement of accuracy.

Fig. 2

Annular summation elements for a circular aperture (xy-plane). (a) Lit region. (b) Shadow region.

Fig. 4

Rectangular summation elements for apertures bounded by straight lines (xy-plane).

The above equations are intended for applications where closed solutions of Eq. (1) cannot be found and will be used in future research. In the remainder of the present paper, their validity will be demonstrated by numerical examples involving the Fresnel diffraction integral where Ugeom is the geometrical field in the absence of diffraction, αF is the modification of the field by diffraction, and where These expressions are valid in the paraxial Fresnel approximation for a distant source point P0 and pertain to a point of observation P = (x,0,z) as in Fig. 1 whereas the point Q = (ξ,η,0) is assumed to lie anywhere in the aperture plane z = 0. In Secs. 2 and 3 of the paper, αF will be reduced to a two-dimensional integral for the respective cases of circular apertures and apertures bounded by infinite straight lines (slit and half plane). The results of the numerical integration will be shown and compared to the corresponding exact solutions, which may be found in Ref. [5].

2. Circular Aperture

For a circular aperture of radius a, Fresnel’s integral in Eq. (4a) can be reduced to a single integral by assuming annular area elements dQ which are centered on the projection Q0 of the point of observation onto the aperture plane, as indicated in Figs. 2a and 2b. With O as the aperture center, x = OQ0 > 0 and ξ = OQ > x for points to the right of Q0, the distance QP defined by Eq. (4b) will then be constant and equal to everywhere on dQ and the integration can be carried out over ξ alone. As the annular area elements are eccentric to the aperture they are, in general, partially obstructed by the aperture rim so that their effective areas will be given by where χ = ∠AQ0Q is the semi-angle subtended at Q0 by the intersection of the area elements with the aperture rim, as indicated in the Figs. 2a and 2b.

When Q0 lies inside the aperture (x ≤ a, as in Fig. 2a), the innermost area elements with radii ξ−x≤a−x are unobstructed and fully contained in the aperture (χ = 0), and the outermost elements with radii ξ−x > a + x are fully obstructed (χ = π). For the intermediate elements the angle χ is found by applying the cosine theorem to the triangle OAQ0 shown in the figure, so that Thus, upon substitution of Eqs. (5a) and (5b) into Eq. (4a) and noting that 2π/(λz) = k/z = u/a2, where the first integral was reduced to an elementary expression by substituting ik(ξ−x)2/2z as the integration variable. When Q0 lies outside the aperture (x≥a, as in Fig. 2b), the inner elements with radii ξ−x≤x−a and the outer elements with radii ξ−x≤x+a are all fully obstructed (χ = π). In the intermediate region, Eq. (5c) applies once again2 and we have

The integrals in Eqs. (5d) and (5e) can now be identified with Eq. (1), with so that, according to Eqs. (3a–c), where The use of these equations on a personal computer is simple. As an example, Fig. 3 compares the numerical values of |α|2 obtained for u = 5π and N = 16 to the exact results given by the Fresnel-Lommel theory. The agreement is good and improves when larger values of N are used, as indicated in the left-hand column of Table 1. The values listed in the table are the maximum errors encountered in the range m≤3N. In this particular case they occurred near the center of the profile (m≤0.5N).

Fig. 3

Approximate (—) and exact (----) diffraction profiles of a circular aperture.

Table 1

The largest errors in diffraction profiles computed from Eqs. (6a), (6b), (6c), (10a), (10b), and (11b).

N	Circular aperture(u=5π)	Slit(u=5π)	M	Half plane
16	6.0×10−2	2.3×10−2	8	1.9×10−2
32	2.3×10−2	5.6×10−3	16	4.3×10−3
64	6.4×10−3	1.4×10−3	32	1.1×10−3
128	2.5×10−3	3.6×10−4	64	2.7×10−4
256	9.6×10−4	9.0×10−5	128	6.6×10−5
512	3.5×10−4	2.3×10−5	256	1.7×10−5

3. Apertures Bounded by Infinite Straight Lines

For a plane aperture of width (l + r) which is bounded by infinite straight lines as in Fig. 4, the reduction of the integral of Eq. (4a) to two dimensions is readily achieved by choosing cartesian coordinates so that the y-axis is parallel to the aperture edges. Letting P = (x,0,z) as before, this leads to where and F(∞) = eiπ/4 is the complex Fresnel integral at infinity.

The amplitude term A(ξ,x,z) in Eq. (1) is now given by the factor that appears outside the integral of Eq. (7b) so that, on letting l = Lh and r = Rh as indicated in Fig. 4 and using Eqs. (3a) to (3d), we find It follows immediately that, if α is known and P is moved to the right or left by ±h, the new value of α is given by the recursion formula which illustrates in a very instructive manner how the diffraction pattern changes when the point of observation is moved, so that one portion of the wavefront is uncovered by the aperture and another portion is covered. The recursion formula Eq. (8b) is convenient for practical computations as it allows the computation of successive values of α without performing the summation of Eq. (8a) for every point. In the examples given below, a starting value for α is obtained from known closed solutions, so that there is no need at all to perform a summation. This use of an exact starting value also improves the accuracy of the computations because it forces the initial error to be zero.

3.1 Slit

For a slit of width 2a, we define l = r = a, L = R = N and h = a/N, so that where it is noted that it is sufficient to perform the computations for m > 0 as the diffraction pattern is symmetrical.

On account of the trigonometric identity this can be transformed into the following expressions, which are convenient for practical computations: As an example, Fig. 5 shows the approximate and exact diffraction profiles of a slit for u = 5π, N = 16, and using the well-known Fresnel solution, , as the starting value. The two curves resemble each other closely, the largest errors being on the order of 0.023 and occurring near m/N = 0.5. The improvement of accuracy achieved by using larger values of N is shown in the center portion of Table 1.

Fig. 5

Approximate (—) and exact (----) diffraction profiles of a slit.

3.2 Half Plane

On letting L = 0 and R = ∞, Eqs. (8a) and (8b) apply to a diffracting straight edge which coincides with the y-axis of Fig. 4. As there is no aperture edge on the right, the last term of Eq. (8b) is now absent and one obtains The previous definition of h as a given fraction of aperture width is no longer applicable but can be replaced by taking h as a certain fraction, say , of the width of the first Fresnel zone. Hence it follows easily that

Figure 6 shows the diffraction pattern computed in this manner for M = 4, using the well-known Fresnel solution, αF(0) = ½, as the starting value. The agreement with the exact result is within ±0.025 for − 4 ≤ m/M < 1, but considerably worse beyond these limits. These discrepancies are reduced when larger values of M are used, as may be seen from the right-hand side of Table 1.

Fig. 6

Approximate (—) and exact (----) diffraction profiles of a half plane.