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.

1. Introduction

The basic theory of Fresnel diffraction at plane apertures was developed long ago [1,2] and is summarized in textbooks [3-6]. For apertures bounded by straight lines (rectangle, slit, half plane), the standard textbook solution in terms of complex Fresnel integrals1 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 semi-quantitative 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.

2. 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 integrals2, or, alternatively, by the Kirchhoff integral,

Here, as indicated in Figs. 1, 2, and 4, 0 is the location of a point source emitting a monochromatic spherical wave of amplitude Asph, circular wave number k = 2π/λ Q is a point in the aperture, dQ is the surface element at Q, is the aperture normal pointing away from the source, and P is the point of observation.

Fig. 1

Notation for rectangular apertures.

Fig. 2

Notation for slits.

Fig. 4

Notation for circular apertures.

In the Fresnel approximation the points P0 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 P0Q 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 where is a small quantity that can be expressed in approximate form. It is well known that, with 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 0 and , and are the distances P0O and OP. The corresponding value of the normal derivative3 in Eq. (2b) is

We now have

The corresponding forms of and UK(P) are essentially the same, except that –cosθ0 is replaced by cosθ and ½ (−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 P0Q 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 is the geometrical field at the point of observation P according to Huygens’ principle, is the modification of the geometrical field by diffraction, E0(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. (2e) are

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 l0, l, … are small compared with unity, so that (l0ξ + m0η)2 and (lξ + mη)2 are much smaller than (ξ2 + η2), and then one finds4

Hence, the magnitude of ϵ (ξ,η) relative to the quadratic term of Eq. (2e) may be estimated as where qmax is the maximum value of (e.g., the radius of a circular aperture) and 〈r〉 is an average of r0 and r. Accordingly, the relative error in Δ(Q) is inversely proportional to the square of the relative distance 〈r〉/qmax. At a distance of ten aperture dimensions, it is on the order of 1 %.

3. Rectangular Aperture

3.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 P0 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 P0P intersects the aperture plane:

This gives θ 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.

3.1.1 Paraxial Approximation

As mentioned in deriving Eq. (4b) above, the direction cosines l0′ and m0′ will be small if the points P0 and P are close to the z-axis of Fig. 1. To a first-order approximation in θ we have l′02, m′02 « 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.

3.1.2 Coordinate Transformation for Off-Axis Sources

According to textbooks, a rotation of coordinates may be necessary when the direction cosines l0′ and m0′ 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 P0P onto the aperture plane. This gives m0′ = 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 P0P, so that l0′ = m0′ = 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 ′ = ′ × ′. Thus, where ϕM and θM are the longitudes and colatitudes of P0 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 θ Then, and where s+ and s− are the same as in Eq. (7b), above. It should be noted that this result differs from the paraxial approximation only by the factors 1/W and W cosθ in the arguments of the Fresnel integrals. Within the above approximation for η′ these factors are equal to unity, and thus Eq. (10b) appears to be no improvement over the paraxial approximation Eq. (7b) of Sec. 3.1.1. It follows that, for rectangular apertures, the coordinate transformations recommended in Refs. [3] and [5] are superfluous. To higher than first order in θ, the ξ- and η-integrals remain inseparable and a closed solution for αF(P) is not possible.

3.1.3 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, y, ϕ and θ 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], where |ϵ (s)| ≤ 2 × 10−3. Accordingly, the following simple algorithm may be used:

Define s.

Let s′ = |s|.

Calculate f(s′), g(s′) from Eq. (11c,d).

Calculate C(s′), S(s′) from Eq. (11a,b).

If s < 0 let C(s) = − C(s′), S(s) = − S(s′). Else, let C(s) = C(s′), S(s) = S(s′).

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].

3.2 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 P0 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 F(t±) = ± ½ (1 + i), [F(t±) − F(t−)] = 1 + i, and Eq. (10b) is reduced to with s as defined by Eq. (7b) but assuming y0 = y = yM = 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 it will also suffice to compute it for positive value values of x, only. The computation is typically carried to a maximum value of θ 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:

Define a maximum (x)max and a step size Δx for x.

Let x = 0.

Compute θ and x from Eq. (12b) and s± from Eq. (7b). Use the algorithm of Sec. 3.1 to find C(s±) and S(s±). Compute αF(P) from Eq. (12a).

Let x = x + Δx. If x < (x)max, go to Step 3. Else, stop.

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. 3

Relative irradiance, |αF(u, v)|2 vs v/u, for a slit. w = 1 mm, z0 = − 2.507 m, z = 1.140 m, λ = 639 nm.

4. Circular Aperture (Fig. 4)

4.1 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(l0, m0, n0) be the geometrical image of P0 at the distance r from the aperture center, so that the position of P relative to C will be given by the vector = r(l − l0, m − m0, n − n0). Let F be the foot of the perpendicular from C onto the xy-plane, let c = FP, and define where

Accordingly, the linear term of Eq. (2e) can be expressed in the form

In the quadratic term of Eq. (2e), we have and, likewise,

In the following, it will be assumed that the points P0 and P are close to the z-axis so that sin2θ0 and sin2θ are negligibly small compared to unity. In this paraxial approximation one obtains

On substitution of Eqs. (14d) and (15c) into Eq. (2e) we have and hence Eq. (3c) is reduced to where the integral over χ was evaluated as 2πJ0(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 (r0 + r)/r0, 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 r0 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, r0 = r = 100 mm are typical and in this case one has u = 0.8p. On the other hand, for limiting apertures used in a radiometer, parameters such as λ = 500 nm, a = 5 mm, r0 = 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).

Fig. 6

Relative irradiances for circular apertures: a, b, c) |αL(u, v)|2 vs v/u for u = 1, 10, and 100 according to Sec. 4.2. d) |αSchw(u, v)|2 vs v/u u for u = 1000 according to Sec. 4.3.

4.2 Lommel’s Solution

Lommel [2] evaluated the integral Eq. (17b) in the form8 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 are Lommel functions of two variables, J(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:

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 so that Airy’s formula, U(P) ∝ J1(v)/v, is obtained from Eqs. (17a) and (3a,b).

For v = 0, we have

For v = u, the well-known relations [9] may be used to show that Accordingly,

The use of Lommel’s equations for numerical computations is straightforward, provided that accurate values of the Bessel functions J(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) or (u/v), provided of course that L and M are evaluated in terms of the Lommel functions V0 and V1 when v < u and in terms of U1 and U2 when v > u.

When v/u or u/v are close to unity, Eqs. (20a–d) will converge on account of the relation J(v) → 0 when n → ∞. The manner in which this limit is approached is illustrated in Fig. 5. As n increases, J(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 V0, V1 or U1,U2, although for consistency it is better to use the former method when v < u and the latter method when v > u.

Fig. 5

Dependence of Bessel functions J(v) on n.

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, are satisfied.

The numerical results presented in this paper were obtained on a personal computer, using standard spreadsheet software9 (a 133 Mhz Pentium computer and Microsoft Excel 7.0). It was found that this software provides accurate values of the Bessel functions J(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 J0(v) from above, as mentioned by Shirley and Datla [13]. Under these conditions, the fine structure of the diffraction pattern poses no difficulties.

The algorithm used in this work for u ≤ 300 was as follows.

Define the value of u, a maximum value vmax, a step size Δv, and the desired decimal accuracy 10−.

Let v = 0. Compute αL(u, 0) from Eq. (22a).

Let v = v + Δv. Compute (u/2) V0(u, v), (u/2)·V1(u, v) from Eq. (20a,b), terminating when Eq. (24a or b) are satisfied. Compute L(u, v), M(u, v) from Eq. (19a,b) and αL(u, v) from Eq. (18b).

If v < u, go to Step 3.

Compute αL(u, u) from Eq. (23e).

Let v = v + Δv. Compute (u/2) U1(u, v), (u/2) U2(u, v) from Eq. (20c,d), terminating when Eq. (24a or b) are satisfied. Compute L(u, v), M(u, v) from Eq. (19a,b) and αL(u, v) from Eq. (18b).

If v < vmax, go to Step 6. Else, stop.

The relative irradiances |αL(u, v)|2 computed in this manner for u = 1, 10, and 100 are plotted in Figs. 6a–c.

4.3 Schwarzschilds’s Approximation

The above-mentioned computational problems encountered with Lommel’s solution when u is large can be avoided by using an asymptotic approximation for αF(P) derived by Schwarzschild [8] in a paper on diffraction effects in defocused telescopes. As this paper is no longer readily available, its contents will be outlined here.

Schwarzschild considered the integral so that his approximation of Eq. (17b) is given by10

The integral W1 is readily shown to be equal to and by evaluating the x -integral of W2 as in Eq. (16b) and then substituting the asymptotic expression [9] one obtains

Letting the first integral in Eq. (26c) can be expressed in the form ∫f dg so that

The second integral in Eq. (26b) can be written as where the first term can again be evaluated by partial integration, using and the second term is a complex Fresnel integral [Eq. (8)]. In this manner, Schwarzschild found

He noted that, by further partial integrations, the resulting expression for W2 would become an asymptotic expansion in negative powers of u but that there was no point in taking the trouble as the last terms in Eqs. (27b) and (28c) are already negligibly small for practical purposes. Therefore,

Schwarzschild estimated that this expression is accurate to 0.005 if u = 100, v/u > 0.2; or u = 300, v/u > 30.

When put in the form of Eq. (25b), Schwarzschild’s approximation becomes 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 = vmin. The choice of vmin can be based on the following table, which was obtained by computing the residuals Δ = |αSchw|2 − |αL|2 for selected values of u.

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.

Define u, vmin, vmax, and Δv.

Let v = 0. Compute αL(0, v) from Eq, (22a).

Let v = v + Δv. Compute (u/2)·V0(u, v), (u/2)·V1(u, v) from Eqs. (20a,b), terminating after the third terms. Compute L(u, v), M(u, v) from Eqs. (19a,b) and αL(u, v) from Eq. (18b).

If v < vmin, go to Step 3.

Let v = v + Δv. Compute δ, β±, s from Eq. (30d). Use the algorithm of Sec. 3.1 to find C(s), S(s). Compute αSchw(u, v) from Eqs. (30b,c).

If v < vmax, go to Step 5. Else, stop.

The values of |αSchw(u, v)|2 computed by this algorithm for u = 1000 are shown in Fig. 6d. They were found to be in excellent agreement with the values of |αL(u, v)|2 obtained from Lommel’s solution.11

It should be emphasized that Schwarzschild’s approximation is different from a superficially similar but significantly less accurate asymptotic approximation of the diffraction integral (17b) cited by Focke [7] and used by Blevin [14], Steel, De, and Bell [16], and Boivin [16] in their work on diffraction errors in radiometry. The respective accuracies of the Schwarzschild and Focke approximations can be assessed by comparing the relative irradiances at the shadow boundary,12 to the exact value |αL(u, u)|2 given by Eq. (23f). From the ratios plotted in Fig. 7 it follows that the Schwarzschild values are accurate to 0.1 % for u > 10, while for u = 100 the Focke values are still off by 6 %.

Fig. 7

Irradiance ratios, |αSchw(u, u)|2/|αL(u, u)|2 and |αFocke(u, u)|2/|αL(u, u)|2, vs u.