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.


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 e ikg() , where k = 2/ is the circular wavenumber of monochromatic light with wavelength .
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, A 1 Ϫ A 2 + A 3 Ϫ ... Ϫ A n , 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 P 0 Q + QP , where P 0 = (x 0 , z 0 ) 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 ( Ϫ x 0 ) 2 << z 0 2 , ( Ϫ x ) 2 << z 2 , and terms in h 2 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 n , as illustrated in Figs. 2a, b, and 4, and define where O is the coordinate origin, Q 0 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 (| z 0 |>> z ) so that 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 1 The value pertains to a point of observation located one aperture halfwidth beyond the shadow boundary. A larger value should be used if necessary.
where u = ka 2 /z is the familiar configuration parameter of Fresnel's diffraction theory for | z 0 |>> 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 mn 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.
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 U geom 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 P 0 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].

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 Q 0 of the point of observation onto the aperture plane, as indicated in Figs. 2a and 2b. With O as the aperture center, x = OQ 0 > 0 and = OQ > x for points to the right of Q 0 , 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 = ЄAQ 0 Q is the semi-angle subtended at Q 0 by the intersection of the area elements with the aperture rim, as indicated in the Figs. 2a and 2b. When Q 0 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 OAQ 0 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 /a 2 , Ϫ iu where the first integral was reduced to an elementary expression by substituting ik ( Ϫ x ) 2 /2 z as the integration variable. When Q 0 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 again 2 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), The use of these equations on a personal computer is simple. As an example, Fig. 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 and F (ϱ) = e i/4 is the complex Fresnel integral at infinity. The amplitude term A ( , x , z ) in Eq. (1) is now given by the factor Ϫ ie i/4 /͙ z 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 ␣ m is known and P is moved to the right or left by Ϯh , the new value of ␣ m is given by the recursion formula  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.

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.
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, ␣ F (0) = (1 Ϫ i )F (͙u /), 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.

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 h = ͙ z /M , of the width of the first Fresnel zone. Hence it follows easily that (11b) Figure 6 shows the diffraction pattern computed in this manner for M = 4, using the well-known Fresnel solution, ␣ F (0) = 1 /2, 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. where = Ϫ x , A = Ϫ ie i/4 /͙z , g ( ) = 2 /z . The only stationary point, defined by g'( ) = /z = 0, occurs at 0 = 0 so that, according to the Appendix of Ref. [2], the stationary-phase value of the integral is where g" ( 0 ) = 1/z and = e +i/4 according as g" ( 0 ) is positive or negative. Since 1/z is positive, this gives ␣ F ≈ 1 so that the "diffraction pattern" of a slit or half plane would be identically equal to one everywhere in the plane of observation.