Objectives of the exercise

•  To experiment with Simpson’s rule for numerical integration;

•  To get experience in the use of SciPy functions;

•  To explore diffraction in the Fresnel and Fraunhofer limits;

•  To gain more experience of plotting data.

Notes:

•  For this exercise, you will be required to submit a program and a short report describing your results.

•  We expect that you will need help as you develop your programs.  Please consult the demonstrators as fre- quently as you need during the drop-in sessions. They are there to help.

Problem: Fresnel diffraction from an aperture – Numerical integration using SciPy

In the computing lectures, you have been briefed about some basic types of numerical integration.   Here we will apply numerical integration to an interesting diffraction problem. The general set-up for the diffraction experiment is shown below. Incident waves, travelling parallel to the z-axis, are diffracted by an aperture, and the scattering pattern observed on a screen. In the Fresnel approximation, the screen is near compared with the size of the aperture, and the diffraction pattern varies with the separation, z. Examples of this “near-field” diffraction are shown in the figure for a square aperture.

A general formula for Fresnel diffraction is given by:

E（x, y, z) =  .-o(o) .-o(o) E（x/, y/) 。xp ╱也(i)z(k) ┌（x - x/)2 +（y - y/)2、、ax/ ay/                                    (1)

where the origin (z = o) is taken at the aperture and k = 也l/λ.  The integral is performed over the x / - y/ plane of the aperture; the electric field of the incident light is taken to be zero everywhere except in the aperture, where it is a

constant, i.e.:

E（x/, y/) = ╱o(E)b

x / , y/ in aperture

otherwise

This is most simply achieved by choosing the limits of integration to fit around the aperture and ignoring the constant

amplitude phase factor of eikz /i:

E（x, y, z) =也(k)l(E)z(b) ..Aperlur(。)e(x)p ╱也(i)z(k) ┌（x - x/)2 +（y - y/)2、、ax/ ay/

(2)

E（x, y, z) is interpreted as the electric field of the diffracted light at coordinates（x, y) on a screen distance z from the aperture. The observed diffraction intensity will be proportional to lE2 l :

I（x, y, z) = ∈b cE（x, y, z)E （x, y, z)*                                                                              (3)

where E* is the complex conjugate of E , c is the speed of light and ∈b is the permittivity of free space.

We need to choose a suitable integration method, and we said in lectures that SciPy has many options available. The Simpson’s method is often considered a workhorse method, which is generally useful.  The implementation in SciPy is easy to use but requires a full set of sample points to be calculated in advance.  The general quadrature approach, in which the name of a function is passed to the integrator together with the limits, is more flexible.  Un- fortunately the library of functions available in scipy.integrate derives from an old Fortran 77 library, which was not designed to take complex data types. (This is strange, because Fortran 77 definitely can handle complex numbers). Since neither approach is ideal, we will try them both!

We will explore the solution to  Eq. (2) in stages.   First we will solve the  1-d diffraction problem,  using both scipy.integrate.simpson() and scipy.integrate.quad() .  Then we will solve the full 2-d diffraction problem using scipy.integrate.dblquad() , which will allow us to calculate solutions for different shapes of aperture.

N.B.Incarryingout the followingstages,youshouldwriteasinglePythonprogram,basedaroundasim-

ilarmenutothatusedinpreviousexercises.Writeaseparatepieceofcodeforeachsectionthatrequiresit,

makingsensibleuseoffunctions.

a)  1-d diffraction: We are going to evaluate the 1-dimensional Fresnel integral. The integral needed is shown in

Eq. (4) both in its complex exponential form, and written out as real and imaginary parts.  Start by writing three functions to evaluate the kernels (inside parts) of the three integrals:

Each function should take the following arguments:

•  the coordinate of a point in the aperture, x/

•  the coordinate of a point on the screen, x

•  the wavenumber, k

•  the distance to the screen, z

To create a diffraction pattern, you will need to consider each point on the screen, x, in turn.  For each value of x (use a loop) you will need to perform the integration over all of the x/ values. To do this you must do the following:

•  For Simpson’s method, create an array of x/ values and another array containing the complex values of the kernel of the integral in the LH part of Eq. (4). Pass these arrays to scipy.integrate.simps() .

•  For the quadrature method, separately evaluate the real and imaginary integrals in the RH part of Eq. (4), by passing the names of your functions and the val- ues of x l(/) and x 2(/) to scipy.integrate.quad() .

Test your program by plotting the intensity against the screen coordinate x. To observe far-field diffraction for a single slit, as shown in the figure on the right, use the parameters listed in the table.  Remember, for the Simpson method you need an odd number of points, N, which implies an even number of intervals.

HINT: To generate the plot, you will need to create two arrays, one for the lE（x)l2 values and the other for the x values. Plot your result using matplotlib as in the previous exercise.

b)  When you are happy that both of your 1-d integration methods are working correctly (they should give identical results with the parameters suggested), explore the effect of changing the aperture width (xl(/)  and x 2(/)) and the screen distance, z, on the appearance of the plot.  Describe and discuss this in your report.  Near-field effects will be observed by increasing the size of the aperture and/or reducing z. You may need to change the screen

coordinate range to accommodate this. You should also examine the effect of varying N on the appearance of your plot from Simpson’s method. Examine the effect of very small N values. Generally, you will need larger N when you explore the near-field features, or when you look far from the axis on the screen.  Can you suggest why?

c)  2-d diffraction: Now you are going to evaluate the 2-d Fresnel integral given in Eq. (2).  However, for general (non-reentrant) apertures this can be rewritten as:

E（x, y, z) = 也(k)l(E)z(b) .y左(/)y/幺.x /左（(x)/幺yp ╱ 也(i)z(k) ┌（x - x/)2 +（y - y/)2、、ax/ ay/                                        (5)

Note that the order of integration is important here. We are taking the y/ integral as the outer integral, with the limits on the inner integral, xl(/) and x2(/)  being potentially functions of y/ .

Proceed in a similar way to part (a), but use scipy.integrate.dblquad() in place of scipy.integrate.quad() , and  don’t  use  Simpson’s  method  for  this  part.    You  will  need  to  provide  two  functions  to  evaluate  the real and imaginary parts of the kernel of  Eq.  (5).   The ordering of arguments is important in these func- tions.   The code snippets  below  indicate  how the arguments should  be arranged and shows  how to call

scipy.integrate.dblquad()

for cases where the inner limits are constant or varying.

To test your approach, choose similar parameters to those in the table above, and calculate the diffraction intensity from a square aperture for a grid of（x, y) values on the screen. You will need to store these values in a 2-dimensional array and plot the array as an image. Here is an example of the use of a 2-d array to generate an image and, on the right, a typical far-field diffraction from a square aperture:

The 2-d diffraction calculations can be slow, so avoid taking too many points in the image; a 上oo × 上oo grid is reasonable for a final image, but will take over 30 s to compute, so you may find more helpful to work with 5o ×5o grids while testing. Repeat your calculation for different values of z and aperture size and shape (square or rectangular). Comment on any qualitative differences in your report.