Algorithms for Fresnel Diffraction at Rectangular and Circular Apertures

This paper summarizes the theory of Fresnel diffraction by plane rectangular and circular apertures with a view toward numerical computations. Approximations found in the earlier literature, and now obsolete, have been eliminated and replaced by algorithms suitable for use on a personal computer.


Introduction
The basic theory of Fresnel diffraction at plane apertures was developed long ago [1,2] and is summarized in textbooks [3][4][5][6]. For apertures bounded by straight lines (rectangle, slit, half plane), the standard textbook solution in terms of complex Fresnel integrals 1 is based on special but poorly documented transformations of coordinates. It is shown in this paper that such transformations cannot be performed accurately for apertures irradiated by arbitrarily located point sources. In the past, the numerical evaluation of the complex Fresnel integrals themselves has also been a problem, and thus previous discussions were confined to simplified special cases. Computational details were omitted and semiquantitative methods (Cornu spiral) were used to describe the nature of diffraction at rectangular apertures.
In the case of circular apertures, the rigorous solution involves Lommel functions of two variables, which are defined as series expansions in Bessel functions and previously had to be evaluated by tedious manual calculations or approximations. For the most part, these approaches have been rendered obsolete by modern computer software. However, approximative methods are still useful for work on personal computers which involves large values of the configuration parameter u defined by Eq. (17a) of this paper. It is shown here that a previously used approximation by Focke [7] is inadequate for this purpose on account of its poor accuracy, but that an older approximation by Schwarzschild [8] gives excellent results.
As algorithms for the computation of Fresnel diffraction patterns on a personal computer have not been published, a compilation of such algorithms is presented in this paper. The underlying theory is stated for off-axis source points, so that the results can be applied to extended sources. For both types of aperture, the closed solutions obtained are paraxial approximations.

The Fresnel Diffraction Integral
The scalar wave function U (P ) associated with diffraction at a plane aperture is customarily expressed by one of the Rayleigh-Sommerfeld integrals 2 , or, alternatively, by the Kirchhoff integral, Here, as indicated in Figs. 1, 2, and 4, P 0 is the location of a point source emitting a monochromatic spherical wave of amplitude A sph , circular wave number k = 2/, Q is a point in the aperture, dQ is the surface element at Q , n is the aperture normal pointing away from the source, and P is the point of observation.
In the Fresnel approximation the points P 0 and P are located at finite distances which are large compared to the wavelength of light and the dimensions of the aperture. Therefore, it is assumed that and that the distances P 0 Q and QP and their normal derivatives do not vary appreciably inside the aperture. Thus they can be replaced, except in the rapidly oscillating exponential function in the integrand, by their values at an arbitrarily chosen reference point O inside the aperture. Under these conditions, the first Rayleigh-Sommerfeld integral Eq. (1a) may be written in the form is a small quantity that can be expressed in approximate form. It is well known that, with O = (0, 0, 0), P 0 = (x 0 , y 0 , z 0 ), Q = (, , 0), P = (x, y, z) (2d) expressed in cartesian coordinates, the required approximation for ⌬(Q ) is where ⑀ ( , ) is the residual error when terms of third and higher order in and are neglected. Here, are the first and second direction cosines of the vectors P 0 O and OP , and are the distances P 0 O and OP . The corresponding value of the normal derivative 3 in Eq. (2b) is We now have U p RS (P ) ≈ ikA sph cos 0 2r 0 r e ik(r 0 +r) The corresponding forms of U s RS (P ) and U K (P) are essentially the same, except that Ϫcos 0 is replaced by cos and 1 /2(Ϫcos 0 + cos ), respectively, where is the colatitude of the point of observation P . Because the diffracted light is confined to a narrow angular range about the central direction P 0 Q unless the aperture dimensions are extraordinarily small, these differences may be judged insignificant. As the Rayleigh-Sommerfeld solutions pertain to the respective cases of p -and s-polarization of the incident light, this implies that Fresnel diffraction is independent of polarization. The Kirchhoff solution has no definable meaning as far as polarization is concerned, but turns out to be equivalent to the Rayleigh-Sommerfeld solutions in the Fresnel approximation. A further solution, the Maggi-Rubinowicz transformation of Kirchhoff's integral [3,4], is not suitable for computations of Fresnel diffraction patterns because it is singular at the boundary of the geometrical shadow.
In this paper, Eq. (2i) will be regarded as the basic form of the Fresnel diffraction integral and will be written as where U 0 (P ) = A sph e ik(r 0 +r) r 0 + r = ͙E 0 (P )e ik(r 0 +r) (3b) is the geometrical field at the point of observation P according to Huygens' principle, is the modification of the geometrical field by diffraction, E 0 (P ) is the normally incident geometrical irradiance at P , and Ϫcos 0 is the inclination factor according to Lambert's law.
The third-and fourth-order terms neglected in Eq.
This shows that the error introduced by neglecting this term depends in a complicated manner on the geometrical parameters involved. Accordingly, it is difficult to assess its magnitude without considering specific cases. However, a few general comments are in order. Under ordinary circumstances, the direction cosines l 0 , l , ... are small compared with unity, so that (l 0 + m 0 ) 2 and (l + m ) 2 are much smaller than ( 2 + 2 ), and then one finds 4 ⑀ ( , ) ≈ Ϫ r 3 0 + r 3 8r 3 0 r 3 ( 2 + 2 ) 2 . (4b) Hence, the magnitude of ⑀ ( , ) relative to the quadratic term of Eq. (2e) may be estimated as where q max is the maximum value of ͙ 2 + 2 (e.g., the radius of a circular aperture) and ͗r ͘ is an average of r 0 and r . Accordingly, the relative error in ⌬(Q ) is inversely proportional to the square of the relative distance ͗r ͘/q max . At a distance of ten aperture dimensions, it is on the order of 1 %.

General Theory (Fig. 1)
When applying the above equations to a rectangular aperture of width 2w and height 2h , it is customary to transform the global cartesian coordinates (x , y , z ) assumed in Sec. 2 into local coordinates (x' , y' , z' ) which depend on the locations of the points P 0 and P and are chosen so that Eq. (3c) is separated into a product of independent Fresnel integrals in and . That is, The first step in this transformation is to place the origin of the local coordinates at the point M where the straight line P 0 P intersects the aperture plane: This gives M being the angle indicated in Fig. 1. The linear term of Eq. (2e) now vanishes and we have This result can be used in two ways to derive a final result.

Paraxial Approximation
As mentioned in deriving Eq. (4b) above, the direction cosines l 0 ' and m 0 ' will be small if the points P 0 and P are close to the z -axis of Fig. 1. To a first-order approximation in M we have l' 0 2 , m' 0 2 « 1, so that the third term of Eq. (6c) can be omitted and Eq. (3c) leads directly to where and is the complex Fresnel integral.

Coordinate Transformation for Off-Axis Sources
According to textbooks, a rotation of coordinates may be necessary when the direction cosines l 0 ' and m 0 ' are too large to justify the paraxial approximation of Sec. 3.1.1, a rotation of coordinates may be necessary. The usual recommendation [3,5] is to place the new x' -axis along the projection of the line P 0 P onto the aperture plane. This gives m 0 ' = 0, so that ␣ F (P ) does indeed assume the form stipulated by Eq. (5). However, the x' -and y' -axes so defined are not parallel to the edges of the aperture, and consequently the two integrals are not separable because the limits of the -integral depend on , and vice versa.
To overcome this difficulty, a different transformation is attempted in the following: The z' -axis is placed in the direction of the unit vector along the line P 0 P , so that l 0 ' = m 0 ' = 0 and Eq. (5) is again satisfied. The x' -axis is chosen so that its projection onto the aperture plane is parallel to the -direction. The y' -axis is defined in the usual manner as y' = z' ϫ x' . Thus, where M and M are the longitudes and colatitudes of P 0 and P with respect to M, and Accordingly, for z = 0, Equation (9d) shows that ' is still not independent of . This was to be expected as it is not possible to rotate the z -axis and have orthogonal x -and y-axes which are both aligned with the aperture edges. Accordingly, the separation of integration limits is not complete unless the first term in the above expression for ' is omitted-a first-order approximation in M . Then, To higher than first order in M , the -and -integrals remain inseparable and a closed solution for ␣ F (P ) is not possible.

Evaluation of Fresnel Cosine and Sine Integrals
The use of Eqs. (7a) and (10b) is straightforward. An example is given in Sec. 3.2, below. It should be remembered that the variation of ␣ F (P ) with P is implicit in Eqs. (7b) and (10b), in that x M , y M , M and M depend on the location of P . It should also be borne in mind that the point M of Fig. 1 will be outside the aperture when P lies in the geometric shadow. This can lead to values of and larger than assumed in Eq. (3a). For this reason, the computation of ␣ F (P ) based on Eq. (7b) or (10b) must not be carried too far into the shadow region.
The only problem that may be encountered on a personal computer is that the Fresnel cosine and sine integrals defined by Eq. (8) are not usually included in standard software packages. For modest accuracy requirements, they can be computed from the equations quoted in Ref. [9], f (s ) = 1 + 0.926 s 2 + 1.792 s + 3.104 s 2 + ⑀ (s ), s Ն 0, (11c) g (s ) = 1 2 + 4.142 s + 3.492 s 2 + 6.67 where |⑀ (s )| Յ 2 ϫ 10 Ϫ3 . Accordingly, the following simple algorithm may be used: If better accuracy is desired, this algorithm can be improved by using the method described in Ref. [10]. Alternatively, software for computing C (s ) and S(s) in Fortran or C can be downloaded [11,12].

Application to Slits (Fig. 2)
The rectangular aperture discussed so far is transformed into a slit of width 2w on setting h = ϱ in Eq. (7b). 5 It may then also be assumed that the source is a long luminous line which is parallel to the slit and passes through the point P 0 in Fig. 1, so that it will suffice to compute the diffraction pattern in the xz -plane shown in Fig. 2. With these assumptions we have t Ϯ = Ϯ ϱ, so that with s as defined by Eq. (7b) but assuming y 0 = y = y M = 0 so that Eq. (9b) is simplified to As the diffraction pattern is centered at and symmetrical about the geometrical source image C shown in Fig. 2, where Fig. 2. Notation for slits. 5 This assumption seems to violate the condition that the values of in Eq. (2e) must be small, but is justifiable on account of Fresnel's zone construction. Because the field at P is not affected by zones located at large distances from the line P0P in Fig. 1, the -integration can be extended to infinity without introducing an error.
it will also suffice to compute it for positive value values of x M , only. The computation is typically carried to a maximum value of M beyond the shadow boundary S , the latter being given by Accordingly, the following procedure may be used to evaluate the dependence of ␣ F (P ) on x: A typical diffraction pattern computed in this manner is shown in Fig. 3. The numerical parameters chosen for this particular example are listed in the figure caption and were taken from an experiment described by Fresnel. 6  (Fig. 4)

General Theory
In evaluating the Fresnel diffraction integral Eqs. (3a-c) for a circular aperture with diameter 2a it is convenient to use spherical coordinates centered at the aperture center O, so that In these coordinates, the path difference ⌬(Q ) defined by Eq. (2e) can be evaluated as follows. Let C = r (l 0 , m 0 , n 0 ) be the geometrical image of P 0 at the distance r from the aperture center, so that the position of P relative to C will be given by the vector CP = r (l Ϫ l 0 , m Ϫ m 0 , n Ϫ n 0 ). Let F be the foot of the perpendicular from C onto the xy -plane, let c = FP , and define and, likewise, In the following, it will be assumed that the points P 0 and P are close to the z -axis so that sin 2 0 and sin 2 are negligibly small compared to unity. In this paraxial approximation one obtains where the integral over was evaluated as 2J 0 (kqc /r ) [9]. On substituting this becomes As expected, these equations describe a circular diffraction pattern which is fully determined by a radial variable, c or v . The pattern is centered at the geometrical source image C , defined by c = v = 0, and ␣ F (P) is constant on any circle about C . The radius of the geometrically illuminated spot at the distance r from the aperture is a (r 0 + r )/r 0 , so that in the notation of Eq. (17a) the geometrical shadow boundary is defined by v = u . The parameter u , which relates the aperture radius a to the wavelength and the distances r 0 and r , 7 can assume widely different values. For example, in the case of a classroom demonstration of Fresnel diffraction, the parameters = 500 nm, a = 0.1 mm, r 0 = r = 100 mm are typical and in this case one has u = 0.8. On the other hand, for limiting apertures used in a radiometer, parameters such as = 500 nm, a = 5 mm, r 0 = r = 1 m are typical, and then one has u = 200. As will be shown later, the diffraction patterns encountered in these different cases are very different. For u → 0 the diffraction pattern approaches the Fraunhofer limit (Airy function), and for u → ϱ it approaches the limit of geometrical optics (rectangle function). (See Sec. 4.2, Figs. 6a-d). 7 It should be noted that, according to Eq. (13c), the distance r = z / cos depends on the location of the point of observation P so that, strictly speaking, u is not a constant if the diffraction pattern is observed at at fixed distance z from the aperture. This dependence of u on P is considered negligible in the Fresnel approximation.

Lommel's Solution
Lommel [2] evaluated the integral Eq. (17b) in the form 8 so that In this notation, is the relative irradiance of the diffracted light and is the phase difference relative to the geometric field. The functions L(u , v ) and M(u , v ) appearing in these equations are defined by where 8 Elsewhere in the literature, L(u , v ) and M(u , v ) are denoted by C(u , v ) and S(u , v ). This practice is not followed here in order to avoid confusion with the Fresnel integrals C(s ) and S(s ) of Eq. (8).
are Lommel functions of two variables, J n (v ) being a Bessel function of the first kind and order n . For checking the accuracy of numerical results, it is useful to note the values of these expressions in special cases: a. In the limit u → 0, Lommel's equations simplify to the familiar Airy formula for Fraunhofer diffraction at a circular aperture. In this case, the entire diffraction pattern lies in the geometrical shadow and Eqs. (20c,d) are reduced to Accordingly, The use of Lommel's equations for numerical computations is straightforward, provided that accurate values of the Bessel functions J n (v ) required for Eqs. (20a-d) are available and the convergence behavior of these equations is taken into consideration.
When v /u or u /v are small, these expansions will converge on account of the monotonic decrease of (v / u ) n or (u /v ) n , provided of course that L and M are evaluated in terms of the Lommel functions V 0 and V 1 when v < u and in terms of U 1 and U 2 when v > u .
When v /u or u /v are close to unity, Eqs. (20a-d) will converge on account of the relation J n (v ) → 0 when n → ϱ. The manner in which this limit is approached is illustrated in Fig. 5. As n increases, J n (v ) exhibits an oscillatory behavior before vanishing after passing through a pronounced final maximum at or below n = v . In this case, L and M can be evaluated in terms V 0 , V 1 or U 1 ,U 2 , although for consistency it is better to use the former method when v < u and the latter method when v > u .
It follows that, for computations of the Lommel functions to m decimals, the expansions Eqs. (20a-d) can be truncated when either of the conditions, ͩ v u ͪ n or ͩ u y ͪ n < 1 2 10 Ϫm , (24a) are satisfied. The numerical results presented in this paper were obtained on a personal computer, using standard spreadsheet software 9 (a 133 Mhz Pentium computer and Microsoft Excel 7.0). It was found that this software provides accurate values of the Bessel functions J n (v ) needed for Eqs. (20a-d) without problems, but that the large number of them required to satisfy Eq. (24b) impeded the speed of program execution when u is large and v ≈ u . In addition, the capabilities of the personal computer were overtaxed by the fact that the diffraction patterns for large values of u are highly structured (see Figs. 6a-d), so that a very large number of data points had to be computed. For these reasons, Eqs. (18a) to (20d) were used in this work only for u Յ 300 while for larger values of the approximation of Sec. 4.3, below, was used. It should be emphasized that this limitation is unnecessary for larger computers. When sufficient computing power is available, the Lommel functions for large values of u can be evaluated efficiently by iterative use of recurrence relations for Bessel functions, beginning at the required large orders and iterating towards J 0 (v ) from above, as mentioned by Shirley and Datla [13]. Under these conditions, the fine structure of the diffraction pattern poses no difficulties. where The use of these expressions on a computer is simple. The only caveat is that they are not valid for small values of v /u so that Lommel's equations must still be used below a suitably chosen minimum value v = v min . The choice of v min can be based on the following For such small values of v /u , only a few terms of Eqs. (20a,b) are sufficient to obtain ␣ L with comparable accuracy. Thus, the following procedure will provide the entire diffraction pattern. Fig. 7. Irradiance ratios, |␣Schw(u , u )| 2 /|␣L(u , u )| 2 and |␣Focke(u , u )| 2 / |␣L(u , u )| 2 , vs u.