## 1.

## INTRODUCTION

Diffraction is the characteristic of light to propagate around obstacles in a limited way, and is a well-known phenomenon in optics.^{1}^{-}^{2} Diffraction influences the real-life behavior of optical systems, such as microscopes and telescopes, and limits the fundamental resolution achievable in these devices. Diffraction phenomena can be classified into two types, near-field (Fresnel) diffraction and far-field (Fraunhofer) diffraction. Fraunhofer, or far-field diffraction is simpler to handle mathematically, because both the wavefronts incident on the aperture and on the observation screen are planar. In contrast, in the case of near-field (Fresnel) diffraction, where either the aperture-source distance or aperture-screen distance is finite, the wavefronts are not planar, and the solution becomes more complex, and the analytical form of the diffraction pattern cannot be found even in the simplest cases. Therefore, numerical methods must be used. Solutions can be derived in terms of Fresnel-Kirchoff, or Rayleigh-Sommerfeld diffraction integrals.^{2} The large amount of numerical data in the calculation of the integrals in real situations means that the use of computers and simulation techniques can be used effectively. With the use of modern computers and software, two-dimensional fast-Fourier transform (FFT) methods can be used to calculate diffraction integrals.^{3}^{-}^{4} For example, the diffraction and propagation of waves from any aperture have been simulated in computers by Rudolf et al.^{4} using Helmholtz–Kirchoff diffraction integrals.

In previous work,^{5}^{-}^{8} we introduced new technique to simulate the diffraction field of apertures having rectangular shapes. This technique is called the Iterative Fresnel Integral Method (IFIM). It uses virtual displacement of the aperture, and repeated calculations of Fresnel integrals to construct the diffraction pattern from rectangular-shaped apertures. The technique can be applied to many interesting optical problems of importance. Using the technique, implemented in computer codes such as MATLAB, the diffraction field from any rectangular-shaped aperture in any situation can be calculated in a few minutes. An amplitude diffraction grating^{1}^{,}^{9} is essentially an array of long parallel rectangular apertures separated by opaque regions. Both the width of the apertures and the opaque regions are of the order of wavelength of light. This grating is widely used in spectrometers to separate the light into constituent wavelengths. In all these cases, the wavefront incident on the grating and that emitted from the grating are both rendered plane by a collimating mirror or lens. In that case, the analytical treatment of the diffraction grating becomes much simpler, because the Fraunhofer diffraction regime can be used to calculate the diffraction fields. The situations, where this Fraunhofer condition is not applocable are not considered or used, because of difficulty of treating Fresnel diffraction from a grating.

In the present Proceedings, we describe how the IFIM technique can be used and extended to the non-trivial case of amplitude diffraction gratings, and discuss some possible applications in teaching diffraction phenomena in classrooms.

## 2.

## INTRODUCTION TO THE ITERATIVE FRESNEL INTEGRALS METHOD

## 2.1

### Theory of a Fresnel diffraction for *N* rectangular apertures

Let us consider the case of *N* rectangular parallel apertures (which is the basis of an amplitude diffraction grating) being illuminated by a light source with a wavelength *λ* at a distance *p* away. The observation screen is located at a distance *q* away. The width of each aperture is *a* and the opaque regions between apertures has a width *b*, the inter-aperture separation being (*a*+*b*). Each aperture has a height of *c* in the *z*-direction. The two edges of the central aperture of the *N*-aperture system (called aperture **0**) are then located at *y _{0}*= -

*a*/

*2*and

*y*’=

_{0}*a*/

*2*respectively, and the edges of the aperture for the next aperture to the right (called aperture +

**1**) are located at

*y*=

_{1}*a*/

*2*+

*b*and

*y*’=

_{1}*3a*/

*2*+b respectively. The next aperture to the right (called aperture +

**2**) are located at

*y*=

_{2}*3a*/

*2*+

*2b*and

*y*’=

_{2}*5a*/

*2*+

*2*b. Finally the edges of the aperture +

*will be located at*

**n***y*=(

_{n}*2n*-

*1*)

*a*/

*2*+

*nb*and

*y*’=(

_{n}*2n*+

*1*)

*a*/

*2*+

*nb*. Similarly, the edges of the first left aperture (aperture -

**1**) are located at

*y*

_{-1}’= -

*a*/

*2*-

*b*and

*y*

_{-1}= -

*3a*/2-

*b*, and, the edges of the aperture

**-**will be located at

*n**y*

_{-n}’ = -(

*2n*-

*1*)

*a*/

*2*-

*nb*and

*y*

_{-n}= -(

*2n*+

*1*)

*a*/

*2*-

*nb*.

As shown in Fig. 1 (shown for the case of *N*=5), the coordinate systems on the aperture and on the image planes are chosen to be centered on the optical axis passing through the center of the central aperture (aperture **0**) and normal to it, and are denoted by (*y*, *z*) and (*Y*, *Z*) axes, respectively. The Huygens−Fresnel principle is then invoked to calculate the total electric field at any given point of the image plane (*Y*,*Z*) by summing up all the contributions of all the elementary wavelets emitted by different area elements inside each of the rectangular apertures. The net complex electric field at the point P (located at the origin of the image plane or screen) due to waves emitted from the central aperture is given by^{1}^{,}^{9}

where *ε _{0}* is the electric field of the source,

*k*is the wavenumber (=2π/λ

_{0}) and

*K*(

*θ*) is the obliquity factor, which accounts for secondary Huygens wavelets emitted from area elements

*dS*in an oblique direction. By integrating the contributions of the Huygens wavelets over the entire area of the central aperture, it can be shown

^{1}

^{,}

^{9}that the total electric field at

*P*is,

where *E _{u}* is the unobstructed electric field at

*P*(i.e. the electric field that would have existed if the entire aperture were absent), and

*C(u)*and

*S(u)*are the Fresnel cosine and sine integrals, being defined by,

Here *w* represents either of the two dimensionless variables *u* or *v*,

The variables *u* and *v* are proportional to Cartesian coordinates *x* and *y*. The limits *u _{0}* and

*u*’ in Eq. (2) are the values of

_{0}*u*corresponding to two left and right edges of the central aperture, i.e. for

*y*=-

_{0}*a*/

*2*and

*y*’=

_{0}*a*/

*2*, respectively. Similarly, the limits

*v*’ and

*v*are the values of the variable

*v*corresponding to the lower and the upper edges of this aperture i.e. for

*z*=-

*c*and

*z*’=+

*c*respectively.

For the rest of these Proceedings, for simplicity, we consider plane-wave illumination of the aperture. This will effectively move the source *S* to infinity, which can be achieved by placing *S* in the focal point of a positive (convex) lens and allowing the collimated light from the lens to fall on the apertures. For the simple case of plane-wave illumination, *p _{0}* is effectively infinity and, equations (4) can be simplified to,

Upon illumination from the source, the total electric field at the center *P* of the image plane consists of the contributions from all of the (*2n*+*1*) apertures. The electric field contribution from the aperture **1** is given, in analogy to Eq. (2), by

Here, the limits *u _{1}* and

*u*’ are the values of the dimensionless variable

_{1}*u*corresponding to two edges of aperture

**1**, i.e. for

*y*=

_{1}*a*/

*2*+

*b*and

*y*’=

_{1}*3a*/

*2*+

*b*, respectively. Similarly, the limits

*v*and

*v*’ are the values of

*v*for

*z*=-

*c*and

*z*’=+

*c*respectively, as before. The electric field contribution by the aperture

**is given by**

*n*where *u _{n}* and

*u*’ are the values of

_{n}*u*corresponding to two edges of aperture

**, i.e. for**

*n**y*=(

_{n}*2n*-

*1*)

*a*/

*2*+

*nb*and

*y*’=(

_{n}*2n*+

*1*)

*a*/

*2*+

*nb*, respectively. The electric field contributions by the apertures on the other side, (e.g., for the -

**1**and –

**n**apertures denoted by

*E*

_{P-1}and

*E*

_{P-n}respectively) are given by equations similar to Eqs. (6) and (7), where all the limits are correspondingly negative. The net complex electric field at P contributed by the all the

*N*(=

*2n*+

*1*) apertures is given by the arithmetic sum of all complex amplitudes

*E*, ….

_{Pn}*E*,

_{P1}*E*,

_{P0}*E*

_{P-1}, …..

*E*

_{P-n}; i.e. by,

Separating the cosine and sine integrals, and using the summation notation for the Fresnel sines and cosines, we can write the total complex electric field at *P* as,

The summation from –*n* to +*n* for both the Fresnel cosine and sine integrals (for the *u* variable only), in effect, carries out the summation of the complex electric field contributions from aperture –* n* to aperture

**, a total of**

*n**N*=(

*2n*+

*1*) apertures in the system. No such summation is involved for the

*v*variable, since the two edges (upper and lower) of the apertures are involved in this

*z*-direction. The net intensity at

*P*is proportional to (

*EE**), i.e.,

where *I _{0}* denotes the intensity of the unobstructed wave due to

*E*, i.e.

_{u}*I*=

_{0}*E*.

_{u}^{2}To calculate the intensity at any off-axis point *P*' on the image plane, one can fix the SOP line and instead of moving the point *P*, one can move the entire aperture system by small amounts in the *yz* plane, so that the relative positions of the aperture edges in the new position and *P* remains unchanged. For example, to find the intensity a point *P*' 1mm below *P*, one can keep the screen undisturbed and move the apertures 1mm upwards, and find the intensity, at *P*, in this situation instead. The point *P* will now see a new set of values for *z*, and therefore, for *v* in Eq.(8). In principle, the intensity at any point *P*' on the image plane can be found in this way by making suitable translations of the aperture in the *y* and *z* directions, and in Eq.(8), substituting the corresponding values of the edges of each aperture in the *u* direction *u _{n}*′..

*n*′.…

_{0}*u*

_{-n}′ and

*u*..

_{n}*n*..

_{0}*n*-

_{n}, (which indicate the positions of the edges of the aperture in the y-direction as seen from

*P*) and

*v*′ and

*v*, (which indicate the positions of the edges of the apertures in the

*z*-direction). This means that new values of Fresnel cosine and sine integrals need to be calculated. Using this method, the entire intensity distribution in the image plane can be mapped out in principle. From Eqs. (8) or (9), it is clear that the calculation of electric field or intensity at a point

*P*requires the evaluation of

*2N*Fresnel cosine and

*2N*Fresnel sine integrals (corresponding to the

*2N*edges of the

*N*aperture system for the

*u*(

*y*) variable), plus two Fresnel cosine and sine integrals for the

*v*(

*z*) variable. The cosine integrals form the real parts of the electric field, and the sine integrals form the imaginary parts. After calculating the complex electric field, the intensity at

*P*is obtained in Eq. (10) by multiplying the field

*E*by its complex conjugate.

_{N}These equations are the basis of calculation of the complete intensity distribution of the Fresnel diffraction pattern from an *N*-aperture, as will be explained next.

## 2.2

### Simulation strategy and methodology

The flow chart of the algorithm is shown in fig. 2. For calculation of the electric field distribution using Eq. (10) at all the points (pixels) on the image plane, a large number of the Fresnel cosine and sine integrals will need to be evaluated quickly and efficiently. This is accomplished in the MATLAB program using the special functions *mfun* (‘*Fresne*lC', *i*:*s*:*f*) and *mfun* (‘*FresnelS*', *i*:*s*:*f*). For example, the MATLAB statement *R*=*mfun* (‘*FresnelC*', *i*:*s*:*f*) generates an array *R* of the Fresnel cosine integrals with arguments staring from *i* and ending in *f*, with a step size *s*.

As shown in fig. 3 for *N*=3 (three apertures), the aperture was displaced (virtually) by an amount *W*, and the extreme limits of the *2N* edges of the displaced aperture was determined and compared to the corresponding limits for the un-displaced aperture to find the range of *u* and *v* values. In a more general *N*-aperture system, when the aperture is moved by W to the left, the ranges for the *2N* edges of the *N*=(*2n*+*1*) apertures will be established as:

for yn: (2n-1)a/2+nb to W+(2n-1)a/2+nb | for yn′: (2n#x002B; 1)a/2+nb to W+(2n+1)a/2+nb |
---|---|

for y1: (a/2+b) to (W+a/2+b) | for y1′: (3a/2+b) to (W+3a/2+b), |

for y0′: (a/2) to (W+a/2), | for y0: (-a/2) to (W-a/2), |

for y-1′: (-a/2-b) to (W-a/2-b) | for y-1: (-3a/2-b) to (W-3a/2-b), |

for y-n′: -(2n-1)a/2-nb to W-(2n-1)a/2-nb | for y-n: -(2n+1)a/2-nb to W-(2n+1)a/2-nb |

Arrays of Fresnel cosine and sine integrals need to be calculated corresponding to these input ranges by the *mfun* statements, with a given step size s. By using these Fresnel arrays, the electric field or intensity dependence in the *y* (*u*) direction can be calculated (first factor in Eq. 9). Following a similar procedure for the entire aperture, i.e. by moving the aperture by *W* downwards, the *z*′ and *z* ranges are determined: *c* to (*W*+*c*) for *z*′ and –*c* to (*W*-*c*) for *z*. Four arrays of Fresnel integrals are to be evaluated for these two ranges, giving numerical values for the calculation of the second factor in Eq. (9). Further details of the computation method and explanation of the MATLAB program is found in Ref. (8) and Ref. (5).

## 3.

## SIMULATION RESULTS

## 3.1

### General

To use the program, the values of the required inputs, i.e. aperture parameters *a* and *b*, height *c*, image area *W*, step size *s*, the aperture-image distance *q*_{0} (all in mm), the wavelength *l* (in nm), the exposure and the number *N* of the slits (apertures) are input through a GUI (fig. 4a). The exposure factor ex is used to adjust the apparent visual intensity in the generated image. The computer generated Fresnel image due to *N*-slits for wavelength *λ*=500 nm, *q*_{0}=400mm, image area *W*=5mm, aperture width *a*=0.1 mm, *b*=0.1 mm (with aperture separation 0.2 mm), *c*=4mm, *s*=0.01mm and *N*=31 is shown in fig. 4b. Since (*a*+*b*)=0.2mm, this system is equivalent to an amplitude grating *N*=31 with 5 lines/mm. The generated diffraction image resembles that of a grating in the near-field, with Fresnel-like characteristics exhibited in the top and bottom edges. Some interference effects can be seen in the left and right edges of the image.

To find the effect of a change of parameters on the diffracted images, we generated a series of images. In the first series, the aperture separation (*a*+*b*) was made successively smaller, while keeping the aperture width (*a*=0.5mm) constant. A series of diffraction images were simulated as shown in figs 5(a)-(d) for a value of *N*=7. In fig. 5(a), in the region between the apertures, no significant interference could be seen between the light from each of the aperture diffraction patterns. But, in fig. 5(b) and (c), some distinct interference between diffracted light from the individual apertures was observed as the separation of the slits is further reduced. As the separation was decreased to zero [fig. 5(d) for *b*=0], a typical single aperture Fresnel diffraction image was obtained for an aperture size 3.5mm×3mm.

In the next series, for *N*=7, the width *a* of the apertures was made successively narrower while their separation (*a*+*b*) was kept constant at 1.5 mm [fig. 6]. From figs 6(a) and (b), as the apertures become narrower, the light is diffracted over a wider region, and some interference between diffracted light can be detected. In figs. 6(c) and (d), at smaller aperture widths, strong mutual interference between diffracted light has generated fringes which look like typical Young-like fringes. Light from each aperture is seen to be diffracted over a wider portion of the image area. In figs 6(c) and (d), fringe density is the same, since this depends only on aperture separation (constant at 1.5mm).

## 3.2

### Comparison with N-slit far-field diffraction

If the aperture-screen distance *q*_{0} is increased, it is well-known in optics that a gradual transition from the Fresnel (near-field) regime to Fraunhofer regime (far-field) should occur. As a general rule of thumb, if *D* is the lateral dimension of an aperture, then Fresnel diffraction occurs if *D* satisfies the following relation^{1}

Otherwise, Fraunhofer diffraction should be observed. In the next simulations, we use a practical diffraction grating with *a*=0.001 and *b*=0.003 (with aperture separation *a*+*b*=0.004 mm, i.e. 250 lines/mm). We selected the following parameters for an 11-slit system, *λ*=500 nm, *c*=2mm, image area size *W* variable, *s*=0.01 mm or 0.05 mm, as appropriate. We then increased the aperture screen distance *q*_{0} in steps from 1 mm to 100 mm while keeping other parameters constant [fig. 7]. A transition from the Fresnel regime to Fraunhofer regime is clearly observed, being consistent with the predictions of Eq. (11). For example, for a 11-slit system with (*a*+*b*)=0.004mm, if we take the total lateral dimension of the slit system to be *D*=0.04 mm, then *D*^{2}/*λ* is about 3.2 mm. We expect the diffraction to be Fresnel-like if *q*_{0} is less than 3.2 mm. If, *q*_{0}> (D^{2}/*λ*)=3.2mm, and we expect the diffraction to be Fraunhofer-like. In between these extremes, a transition from Fresnel to Fraunhofer regime should occur [as in figs 7(c) and 7(d)]. The usual mathematical analysis of an *N*-slit system^{1,9} deals with the N-slit pattern only in the far-field Fraunhofer regime. The far-field diffraction pattern of an *N*-aperture system in the image plane (*Y*, *Z*) can be exactly calculated by an analytical formula^{1,9}

where *β*, *γ* and *α* are defined as,

Here, *a* is aperture width, (*a*+*b*) is the inter-aperture separation, *c* is aperture height, *λ* is the wavelength and *R* is the aperture-screen distance (assumed to be sufficiently large). As stated above, in the computer simulations, this far-field conditions should be amply satisfied in figs 7(e) and 7(f), for *R* (or *q*_{0})=50 mm and *R*=100 mm, respectively.

In the usual analysis of far-filed Fraunhofer diffraction from an *N*-slit system,^{9} it is shown that the *principal maxima* of the diffraction pattern occurs when the [*sin Nγ*/*sin γ*]^{2} factor have maximum values, i.e. when,

Here, *θ* is the diffraction angle and *m* is the diffraction order. For *m*=0, we have the zeroth-order, for *m*=1, the first order, and so on. In addition to these principal maxima, several much weaker *secondary maxima*, and several minima of zero intensity are predicted to occur with equal separation between the principal maxima. If *N* is the number of grating slits (apertures), then the number of secondary maxima is equal to (*N*-2), and the number of minima is equal to (*N*-1).^{9} To make comparisons, we reproduce fig. 7(f) next, with the intensity values multiplied by 10 to bring out the faint details. We observe 7 principal maxima, with *m*=0 (at *Y*=0), with *m*=+1 (at *Y*=12.5 mm), *m*=+2 (at *Y*=25.03 mm), *m*=+3 (at *Y*=37.59 mm), with *m*=-1 (at *Y*=-12.5 mm), *m*=-2 (at *Y*=-25.0 mm), and *m*= -3 (at *Y*=-37.5 mm). In addition, 9 (i.e., *N*-2) secondary maxima can be discerned between the principal maxima at the zeroth order and the first and the second orders, along with 10 minima of zero intensity, interspersed between them. The values of sin *θ* can be estimated by finding the ratios (*Y*/*R*, with *R*=100 mm) from the computer-simulated image (Fig. 8). These calculated values are: 0.125 (for *m*=+1 and *m*= -1), 0.2503 (for *m*=+2 and *m*= -2), and 0.375 (for *m*=+3 and *m*= -3) From Eq. (14), on the other hand, for *λ* =500 nm, (*a*+*b*)= 0.004 mm, and the following values of sin *θ* are calculated: 0.125 (for *m*=+1 and *m*= -1) 0.250 (for *m*=+2 and *m*= -2), and 0.375 (for *m*=+3 and *m*= -3). We therefore conclude: the computer-simulated results of the diffraction field, in the far-filed Fraunhofer limit, precisely agree, both qualitatively and quantitatively, with the results from the exact Fraunhofer theory, as represented by Eqs. (12) and (14). This agreement gives strong support to the validity and correctness of the present simulation methods.

## 4.

## DISCUSSIONS

In this Proceedings, we describe a technique to simulate the near-field diffraction image of a *N*-aperture system or, equivalently, an amplitude diffraction grating, using the iterative Fresnel Integrals method. The theoretical background and the implementation of the algorithm are described. Some examples of computer-simulated images as well as the practical case of a 250 line/mm amplitude diffraction grating are presented in the results. For the simulations, an Intel core-i7 PC having 4GB of RAM and running at 3.4 GHz was used. The simulation time for a typical 2000X2000 pixel image for *N*=11 is less than 45 seconds. The time depends on the computer configuration used, or, on the image size. The expected transition from Fresnel to Fraunhofer regions can be clearly observed under the appropriate conditions. For example, when the aperture-screen distance *q _{0}* is increased sufficiently into the far-field regime, the predicted principal maxima in the Fraunhofer limit are observed at the expected positions for different diffraction orders, and the predicted number of secondary maxima and minima agrees with the far-field Fraunhofer theory. The simulation techniques and the resulting virtual experiments provide a simple and easy-to-use method to study the complex phenomenon of diffraction from amplitude diffraction gratings. It can be used to teach diffraction phenomena from amplitude gratings in a classroom setting, where the virtual experiments can be performed in real-time in front of the students.

## REFERENCES

Hecht, E., [Optics], 4th ed. Pearson Education, Singapore, Chap. 10 (2002).Google Scholar

Born, M. and Wolf, E., [Principles of Optics], 7th ed., Cambridge University Press, Chap. 8 (1999).Google Scholar

Wilson, R.G., McCreary, S. M. and Thompson, J. L., “Optical transforms in three-space: Simulations with a PC,” Am. J. Phys. 60 49–56 (1992).Google Scholar

Rudolf, P.G., Tollett, J.J. and McGowan, M.M., “Computer modeling of wave propagation with a variation of the Helmholtz-Kirchhoff relation,” Appl. Opt. 29, 998 (1990).Google Scholar

Abedin, K.M., Islam, M.R. and Haider, A.F.M.Y., “Computer simulation of Fresnel diffraction from rectangular apertures and obstacles using the Fresnel integrals approach,” Opt. Laser Technol. 39, 237–246 (2007).Google Scholar

Abedin, K.M. and Rahman, S.M.M., “Computer simulation of Fresnel diffraction from double rectangular apertures in one and two dimensions using the iterative Fresnel integrals method,” Opt. Laser Technol. 44, 394–402 (2012).Google Scholar

Abedin, K.M. and Rahman, S.M.M., “The Iterative Fresnel Integrals Method for Fresnel Diffraction from Tilted Rectangular Apertures: Theory and Simulation,” Opt. Laser Technol. 44, 939–947 (2012).Google Scholar

Abedin, K.M. and Rahman, S.M.M., “Fresnel Diffraction from N-apertures: Computer Simulation by Iterative Fresnel Integrals Method,” Optik 126(23), 3743–3751 (2015).Google Scholar

Jenkins, F. A. and White, H. E., [Fundamentals of Optics] 4th ed. McGraw-Hill, New York, 358–362 (1976).Google Scholar