Part II Computaonal Physics Exercises
David Buscher v. – //
Contents
Introducon
Assessment
. Pigeonholes ………………………….. . Codequality…………………………..
Geng started
. UsingtheLinuxMCS……………………… . Sengupyourscipyenvironment ………………. . Sengupyourworkspace…………………… . Runningyourprogram …………………….. . Plongyourresults………………………. . Examplecode ………………………….
Exercise : Integraon
. Goal………………………………. . BackgroundTheory………………………. . Tasks ……………………………… . Stepbystepguide:Monte-Carlointegraon . . . . . . . . . . . . . . Step by step guide: Numerical Integraon with . . . . . . . .
Exercise : The Driven Pendulum
. Goal……………………………….
. Physics…………………………….. . Tasks ……………………………… . Hints ………………………………
Exercise A: Helmholtz Coils
. Goal………………………………. . Physics…………………………….. . Tasks ……………………………… . Hints ………………………………
Exercise B: Diffracon by the FFT
. Goal………………………………. . Physics…………………………….. . Tasks ……………………………… . Hints ………………………………
Introducon
These exercises are designed to help you learn some , understand some physics and learn a lile more about numerical programming. You should aempt the exercises one per week, in the order given. Demonstrators will be on hand in the MCS to assist.
The speed at which people can program varies widely. It even varies for individuals day by day — all programmers will sympathise with how much me can be wasted trying to fix an obscure bug. So do keep an eye on the clock, and do seek assistance. Talk to your colleagues so you can share experience and spot each others’ mistakes. The goal of the exercises is for you to learn – and you can only do this by experience.
Each exercise is split into tasks. The core task or tasks are ones that I hope you will all be able to achieve in the class in the me available. If you hand in a working soluon to this you will get approximately –% of the credit for the exercise. There is also one or more supplementary tasks per exercise, which take the problem further, or in a different direcon. I hope that many of you will achieve some or most of these, and they count for the remaining fracon of the credit.
I have given structured hints to each problem. In addion, you should browse the code examples available on the MCS computers (see §.), and, importantly consult the relevant documentaon for the libraries such as and — finding and learning from online documentaon and code examples is an important skill. There are some starng links from the course web-page or just use a search engine.
Assessment
Over the weeks you should carry out three exercises: Exercise , Exercise , and either Exercise A or B. The marks available are a maximum of for exercises and , and a maximum of for exercise A or B, for a total of marks.
You should upload your soluon to the exercises to your pigeonhole (see below). The deadline is posted on the course web site and the TiS. You are encouraged to submit your work as soon as it is in a reasonable state, and move on. So don’t spend too much me on the tasks, but don’t spend too lile me either.
The exercises suggest what to submit: normally this is
. Source code file: normally this will be just one or two files (or files). Python files should run without problem on the MCS using ). If using , your files should compile and link on the MCS system (I will assume that the files require the Gnu Science Library ( ) and and libraries to link).
. a text file which describes very briefly what you did, any major problems, and menons any numerical answers required and describes any accompanying plots. An alternave is to submit a PDF file, called which contains the text menoned above, but can also include the plots men- oned below. Do not submit Microso Word files or similar — only plain text or PDF files will be considered.
. PDF plots illustrang the problem. See §. on how to plot in PDF format. Your plots can be included in the file if desired, in which case addional plot files are not required. Remember to annotate your graph axes appropriately.
Please be aware of the following:
• You must submit Python files and not files (Jupyter notebook files). The Jupyter notebook is a good environment for experimenng with Python, but for this course, you must submit your work as a standard Python program.
• Please do not submit data files or C executables to the pigeonholes. Any files (with the excepon of PDF files) larger than MB will be deleted without looking at them.
. Pigeonholes
You are required to submit your work electronically to what I have called your pi- geonhole. You should upload your answer to each exercise to its own directory. These directories should have the names
and either or
respecvely, where should be substuted with your own crsid.
Important: Let me know if you have problems using your pigeonhole or if it is missing. If you find you can’t delete files once uploaded, and need to add a new version, just upload a new version with a number at the end, e.g. Or you can create a whole new directory in the pigeonhole and put a whole new soluon there; perhaps call it or something like that.
. Code quality
In the first two exercises, the main emphasis of the marking will be on successful compleon of the tasks. In the last exercise a percentage of the marks will be given for code quality (this will also be true for the oponal project). High-quality code will include at minimum:
• Structured code: tasks broken down into sensible funcons;
• Meaningful funcon names;
• Meaningful variable names (this is less important than for naming funcons:
single-leer variable names can be meaningful if the meaning can be inferred
from context, e.g. loop counters);
• Appropriate levels of commenng (at minimum idenfying what each func-
on does);
• Sensible use of whitespace to indicate code structure.
Remember, code quality is important because (a) it helps to make the code easier to understand and debug for the person wring the code and (b) it makes the code easier to understand and debug for others who may have to modify the program later on.
Wring high-quality code is somewhat like wring high-quality prose or mathe- macs: clearly structuring and explaining your thoughts for someone else can help clarify your own thinking, and this clear explanaon is what code quality is all about.
Geng started
. Using the Linux MCS
If you are not familiar with using the MCS Linux system, I strongly recommend that
you read the relevant secons of the Part IB compung tutorial available at
(this docu- ment is also linked from the Part II Computaonal Physics course website), specif- ically secon . about logging in and using a text editor, and the appendix with a list of useful Unix/Linux commands.
. Seng up your scipy environment
Unfortunately the Python setup on the MCS is not updated very frequently, and parcularly the library is quite out-of-date. The best way to avoid issues arising from this is to install your own copy of scipy, by typing the following into the terminal at the command prompt:
This could result in a fair amount of acvity on the screen before returning you to a command prompt. You can then verify that you have a recent version of scipy by typing the following command into your terminal:
This prints out the version number; if your version number is at least , then
the examples from the lectures should all work. . Seng up your workspace
I strongly recommend using a new directory for each project you do — this will make it much easier to look aer your files. To do this, use the UNIX command. So for example from your home directory:
Then start work in this directory; next week, make etc. . Running your program
You can use a text editor such as , or to type in your code. Once you have saved it in, say, you can run it by typing
andlookingattheresults.Notethatifinsteadyouweretotype , this would invoke the interpreter. In this course we are emphasising the use of , which is very similar to, but not the same as, .
. Plong your results
It’s oen convenient to use to write output from a program. You can catch this in a file using UNIX re-direcon:
You can then read the answers back in to a separate program for plong. Alter- navely, you can skip saving the results to a file and plot the results from the same program which calculates the data. This is ok as long as the calculaon me is short (less than a few seconds): otherwise when you make a mistake in the plong you need to wait for the calculaon to repeat.
Remember to hand in plots in PDF format. Be sure to add capons for the axes and a tle. Here is a simple Python program which plots a sine funcon with axes and saves it as a PDF file.
. Example code
There are example programs which might be useful in the directory:
.
Remember the command which searches through files for text — you
could use this to search files in this directory for example.
Exercise : Integraon
. Goal
To learn how to evaluate integrals numerically, using (a) a self-wrien Monte-Carlo
method and (b) a general purpose integrator from .
. Background Theory
Suppose we have a funcon f of n parameters which we can regard as an n- dimensional vector r. We want to integrate this over a muldimensional volume V. We can esmate this by taking N samples distributed at random points ri through- out the volume V, as follows:
where N−
fdV≈ V⟨f⟩±σ ()
⟨f⟩ ≡ N
⟨f⟩ − ⟨f⟩ /
f(ri) () σ≈V N ()
and
with
i=
N−
⟨f⟩ ≡ N f(ri) ()
i=
Note that the error esmate in equaon () is not guaranteed to be very good, and is not Gaussian-distributed: treat it as indicave only of the error. In this exercise you will esmate the error by a robust Monte-Carlo approach and compare it with the theorecal error.
. Tasks
Core Task : Write a program to find an approximate value of this integral and an associated error esmate:
s s s s s s s s
sin(x+x+x+x+x+x+x+x)dxdxdxdxdxdxdxdx
where
s=π
Show that the error in the integral falls off as N−/ where N is the number of Monte-Carlo samples. In this task you should esmate the error on the value for a given N from the standard deviaon of several independent esmates, and plot a suitable graph.
(In case you are wondering, there is an analyc answer in this case, but that’s not the point! The answer is in fact ×(− sin(π/)+ sin(π/)− sin(π/)) ≈ ..)
using Monte-Carlo techniques.
Upload to the course pigeonhole: Source code ( or ), relevant graph- ical plots, and a containing a couple of sentences summarising what you managed to achieve, and the integral’s value and the error es- mate. Specifically, include your integral value, error esmate and N value for the largest value of N for which you can compute results in a reasonable amount of me.
Note: if you are wring in , then wring a vectorised program will allow you to reach much higher values of N than if you code all the loops explicitly. Try this if you have me.
Core Task : Write a program to evaluate the Fresnel integrals accurately using a standard integraon roune from the , and use this to make a plot of the Cornu spiral using . Do not use a Monte-Carlo roune — they are not effi- cient for low-dimensional integrals — instead use a standard quadrature technique. One version of the Fresnel integrals can be wrien
u πx C(u)=cos dx
u πx S(u)=sin dx
Supplementary Task Compare your error derived from the scaer of your Monte- Carlo simulaons with the theorecal error esmate given by equaon (), using a suitable plot.
Supplementary Task : Use the previous work to calculate and plot the diffrac- on paern of a slit of width d = cm illuminated at normal incidence by plane monochromac waves of wavelength λ = cm. The paern is observed on a plane normal to the slit and the incoming waves. Calculate the paern when the screen is at distances D = , and cm from the slit. Plot both the relave amplitude and the phase of the paern. (The absolute amplitude and phase are not important here).
Physics reminder: the complex amplitude in the near field on the axis of an open slit
Upload to the course pigeonhole: source code, relevant plots.
Upload to the course pigeonhole: source code, relevant plot.
is given by
x
Ψ ∝ (cos(πx /(λD)) + i sin(πx /(λD))) dx x
so you can use the Fresnel integrals directly scaling the length dimensions by /(λD). This result holds good for D ≫ d /, And the Fraunhofer limit applies when
D ≫ d . λ λ
Upload to the course pigeonhole: source code, relevant plot.
. Step by step guide: Monte-Carlo integraon
. The first goal is to write a funcon that returns an esmate of the integral’s value using N sample points using Equaon . It is sensible to make N an argument (perhaps the only argument) of this funcon, and make the funcon return the value of the esmate. This funcon can then be called over and over to get esmates of the integral.
. In , the value of π is held in . You can simplify your life a lile byhavinga statementatthestartofyourcode,then you only need to type thereaer to get its value.
. You will need to be able to generate random points in -dimensional (or gen- erally n−dimensional) space. Use the random number genera- tors for this. Look in the example code from the lectures and in the examples directory to see how this can be done.
. The basic logic you require now is to loop over increasing values of N; and for each N value, esmate the integral several mes, and find a best value and error by looking at the mean and spread of esmates returned. How big should N be? Experiment! And think about the error you might want to achieve. To debug your code, a small value is good, but when you are sure things are working well you should set N to be as large as possible.
How many samples can you compute in a few minutes on the MCS? Does this makes sense, supposing the computer clock speed is about GHz say? In , your code may run orders of magnitude faster if you think about how to vectorise your code, in other words use funcons such as to act on a whole array of numbers in a single funcon call.
. For a given N, esmate the integral nt mes. (nt = might be a sensible choice). (Make sure your random numbers are different for each trial – they will be as long as you do not re-seed the random number generator). From these nt independent esmates you can esmate the best value (from the mean of the values) and the error in the value (from the rms of the values). You could either store these independent samples in an array so that you can compute mean and standard deviaons at the end, or you can compute the mean and mean square values as you go through your loop over nt.
. Make your program output the esmates and errors, and plot them. For the longer runs it is likely best to save the results to a file and read them back in to a separate plong program.
You could perhaps use the funcon in to illustrate
error bars. Also don’t forget the logarithmic plong capabilies in
for example.
. Save your plot to a file for submission to your pigeonhole.
. Step by step guide: Numerical Integraon with
. You will need to choose a roune to do the integraon. The funcon is a good general purpose integrator. Example code is available online: search for “scipy numerical integraon examples”.
. Write suitable funcons that evaluate the two integrands. Make sure the funcons have the correct parameters and return value types to match the funcon required by the relevant integraon funcon.
. Evaluatetheintegralsforvariousvaluesofsanduse toplotthespiral.
. Plong the diffracon paern should now be straighorward.
Exercise : The Driven Pendulum
. Goal
To invesgate the physics of a damped, driven pendulum by accurate integraon of its equaon of moon.
. Physics
The pendulum comprises a bob of mass m on a light rod of length l and swings in a uniform gravitaonal field g. If there is a resisve force equal to αv where the bob speed is v, and a driving sinusoidal couple G at frequency ΩD, we can write
m l dθ = −m g l sin(θ) − α l dθ + G sin(ΩDt) () dt dt
Rearranging:
where we have defined q ≡ α/(m l) and F ≡ G/(m l).
dθ = −g sin(θ) − q dθ + F sin(ΩDt) () dt l dt
For the purposes of this exercise, take l = g, so the natural period for small os- cillaons should be π seconds. Also let the driving angular frequency be Ω = / rad s−. This leaves q and F, and the inial condions, to be varied. For all of these problems, we will start the pendulum from rest i.e. θ ̇ = at t = , but we will vary the inial displacement θ. The parameter space to explore then is in the three values q, F and θ.
. Tasks
Core Task First re-write this second-order differenal equaon Equaon as a pair of linked first-order equaons in the variables y = θ and y = ω = θ ̇ . Now write a program that will integrate this pair of equaons using a suitable algorithm, for example a th order Runge-Kua technique, from a given starng point θ = θ and ω = ω at t = . I recommend using rather than implemenng your own Runge-Kua technique. The program should write out line by line the values of t, θ, ω and total energy at frequent enough intervals to follow the evoluon of the equaons. If your files get big, you should modify the program to write out the results less frequently.
Test the code by seng q = F = and starng from (θ,ω) = (.,.), and plong the soluon for , , …natural periods of oscillaon. Overlay on your plot the expected theorecal result for small-angle oscillaons — make sure they agree!
Test how well your integrator conserves energy: run the code for, say, , os- cillaons and plot the evoluon of energy with me.
Now find how the amplitude of undriven, undamped oscillaons depends affects the period. Plot a graph of the period T versus θ for < θ < π. In your wrien answer file, state the period for θ = π/.
Core task Now turn on some damping, say q = , , , plot some results, and check that the results make sense. Now turn on the driving force, leaving q = . from now on, and invesgate with suitable plots the behaviour for F = ., ., ., .. What happens to the period of oscillaon? Note that the period of oscillaon is best observed in the angular velocity rather than angular posion, to avoid problems with wrap-around at ±π.
Supplementary task Invesgate the sensivity to inial condions: compare two oscillaons, one with F = ., θ = . and one with F = ., θ = .. Integrate for a ’long me’ to see if the soluons diverge or stay the same.
Upload to the course pigeonhole: Source code, containing a couple of sentences summarising what you managed to achieve, and the value of the period for θ = π/; include an appropriately-scaled plot showing how well energy is conserved in at least one case and a plot of period versus amplitude.
Upload to the course pigeonhole: Source code, a sentence or two in the
explaining what you see in the soluons, and illustrave plots
of the displacement θ and the angular velocity dθ versus me. dt
Upload to the course pigeonhole: Source code, a sentence or two in the explaining what you see in the soluons, and illustrave plots of the behaviour.
Supplementary Task Try plong angle versus angular speed for various solu- ons, to compare the type of behaviour in various regimes: you can invesgate chaoc behaviour using this simple code — have a look in the books or a web site for examples. There is a nice demo for example at . There’s lots more physics to be explored here — experiment if you have me!
. Hints
. Hopefully rewring the equaon as ODEs is straighorward: if not, ask!
. I recommend that you use the fourth-order Runge-Kua integrator for ex- ample the funcon, rather than wring your
own.
. Both techniques are best done by first wring a funcon that evaluates the
derivaves of the ODEs at a given me given the current values of the vari-
ables. You can look at the example programs in or which solve the spinning ring problem discussed in lectures. You can perhaps
adapt some of that code if you get stuck.
. Use the ”>” redirect symbol on the command line to send your output to a file. It is probably not necessary to write out every me step used by the integrator (this is not a problem with the integrator as it only returns values at the ordinates you choose): you may want only to print out every, say, th or th value.
. A note on plong: you may want to write out a text file containing columns of data — perhaps me, angle, angular speed, energy. You can then read back the results into a separate Python program and plot them over different me spans.
. To find the period versus amplitude relaonship, a simple (and just about acceptable) way is to measure the period off a suitable plot, and do this for several values of θ. However, it is much beer to alter your code to esmate the period directly. You can then loop over θ values and plot the period versus amplitude relaon. Two obvious approaches spring to mind. First, you can find when θ first goes negave. This is when the me is approximately T/ where T is the period. How accurate would this result be? You could also find several zero crossings by considering when y changes sign or becomes
Upload to the course pigeonhole: Source code, a sentence or two in the explaining what you see in the soluons, and illustrave plots of the behaviour.
exactly zero aer a step is taken; by counng many such zero crossings and
recording the me between them you can get a more accurate value for T. . Note that you need to think about what happens when the pendulum goes “over the top” and comes down the other side — you need to think about the π ambiguies involved, and what this means in terms of the “period” of an
oscillaon.
Exercise A: Helmholtz Coils
. Goal
Write a program to calculate the magnec field caused by Helmholtz coils. Check your soluon agrees with the on-axis analycal result. Invesgate with suitable plots the uniformity of the field near the centre of the system.
. Physics
The physics of this exercise is straighorward: you will all know that a small length
of wire dl carrying a current I creates a magnec field
dB= μ Idl∧r ()
π r
at a locaon r with respect to the current element. So by breaking up any wire into short elements we can simply add up all the contribuons to find the total magnec field. Note that in this problems you are expected to do simply this. There are of course more accurate ways to compute integrals numerically — but this exercise is more about how to organise your code than how to evaluate integrals to high accuracy.
Recall that the axial field on the axis of a single coil is given by
B = μ I R ()
(R + x)/
Helmholtz coils produce a fairly uniform field close to the centre of the system. They consist of two co-axis coils carrying parallel and equal currents I. The separa- on D of the coils centres is set equal to their radius R. For coils with axes aligned with the x axis we can place the coil centres at (±(R/), , ). For the purposes of this problem set R = . m and it is convenient to set the current so it has a (very large!) value of (/μ) amps, to get easy to interpret numbers (although somewhat high field strengths!). The field expected at the coil centres is then
B = / μ I R
. Tasks
Core Task Write a program which computes the magnec field from a single coil of radius R = m carrying a current I = /μ amps. Put the coil centre at (, ) and let the coil’s axis run in the x direcon. Calculate the field in the x − y plane. First check your field on axis (y = z = ) agrees with theory by plong your result against the theorecal one. Then calculate the field on a x − y grid and plot your results. To compare theory with your calculaon, you should plot the difference between theory and calculaon along the x axis, as well as the actual values.
You will need to break the wire up into a series of straight line elements, and using the Biot-Savart law, you can then add up the total field. How many elements should you break your circle into? You should experiment with different values in your code. And how does the final accuracy of your calculaon depend on this value?
Core Task Now adapt your program to calculate the field from Helmholtz coils arranged as in the diagram. Plot the field strength near the centre. How uniform is the field? The coil centres should be set to be (±. m, , ) with the coils’ axes poinng in the x direcon. Calculate and plot the magnec field strength in the vicinity of (, , ). Show that the field is quite uniform in this region. To quanfy this uniformity, find the maximum percentage deviaon of the field magnitude from that at (, , ) within a cylinder which has diameter cm, length cm, and is coaxial with the coils. The centre of the cylinder should be at (, , ). Because of the symmetry, it is of course sufficient to calculate and plot the field for z = i.e. plot B(x, y, ), with x and y extending to ±cm.
Supplementary Task Extend or modify your code to invesgate the field caused by a series of N coaxial coils with uniform spacing carrying the same current. Make plots to show the effects of varying N while keeping the distance between the outermost coils constant at D = R.
. Hints
There are several ways to approach this problem. Do some thinking before you start coding. First, you have to decide how to deal with the vector fields, both for the posion and magnec field vectors. One possibility is to have three arrays represenng the Cartesian x, y and z components of the vectors respecvely. Al- ternavely, a higher-dimensional array where the last dimension of the array is of length three may be the way forward. The choice is yours.
You then need to think about how to represent the wire. I recommend that you write a program that is general, i.e. that can calculate the field due to an arbitrary wire. You could define the wire by a set of N line segments defined by N + points in space, each secon of wire running from (xi, yi, zi) to (xi+, yi+, zi+), i.e. from ri to ri+ using vectors.
You can now write a funcon that calculates the contribuon dB from each secon using equaon Equaon , and add them up to find the total field of the wire. You need to evaluate the cross product of course. Be careful close to the wires where r gets small.
Now you need to create a wire which represents a single coil suitably, and evaluate the field it produces by stepping over a -d grid of points in the (x, y) plane.
You can plot the -dimension variaon of the magnitude of the field as an image using the funcon from . You could also experiment with plong the vector field using funcon. There can be a lot of data, so it may be simplest in some cases to compute your results and plot them in the same program, avoiding the complexity of wring out and then reading back in the data.
Exercise B: Diffracon by the FFT
. Goal
Write a program to calculate the near and far-field diffracon paerns of an ar- bitrary one-dimensional complex aperture using the Fast Fourier Transform tech- nique. Test this program by using simple test apertures (a slit) for which the theoret- ical paern is known. Invesgate more complicated apertures for which analycal results are difficult to compute.
. Physics
Plane monochromac waves, of wavelength λ, arrive at normal incidence on an
aperture, which has a complex transmiance A(x). The wave is diffracted, and the
Figure : Geometry for diffracon calculaon
paern is observed on a screen a distance D from the aperture and parallel to it. We want to calculate the paern when the screen is in the far-field of the aperture (Fraunhofer diffracon) and also in the near-field (Fresnel).
Using Huygen’s construcon, we can write the disturbance at a point P on the screen, a distance y from the axis, as
ψ(y) ∝ A(x) exp(ikr) dx r
where k = π/λ. We have assumed that all angles are small: x, y ≪ D
so that we are close to the straight-through axis and can therefore neglect terms like cos(θ) which appear if we are off-axis. We now expand the path length r in
powers of x/r:
r =D +(y−x)
y xy x (y−x)
r≈D+D− D +D+O D
If we now neglect the variaon in r in the denominator of the integral, seng r ≈ D,
which is adequate for x, y ≪ D, then we can write
ψ(y)∝exp(ikD)expiky A(x)expikxexp−ikxydx ()
D D D D
The diffracon paern is thus the Fourier-transform of the modified aperture func-
on A′: with
ψ(y)∝expiky A′(x) exp−ikxydx () D D
A′(x) = exp ikx A(x) () D
In the far-field (Franhofer limit) we have kx/(D) ≪ π so that A′ ≈ A for all values of x in the aperture where A(x) is non-zero, i.e. the familiar result
d ≫ xmax . λ
The distance xmax/λ is the Fresnel distance. In this case, the diffracon paern is just the Fourier transform of the aperture funcon.
Note that we can calculate the near-field (Fresnel) paern also if we include a step to modify the aperture funcon according to Equaon before we take its Fourier transform.
Note that if we are only interested in the paern’s intensity, we can ignore the phase prefactor in Equaon .
Finally, we can discreze Equaon , by sampling the aperture evenly at posions
xj
where Δ is the distance between the aperture sample posions xj.
N− −ikxj y
ψ(y) ∝ Δ A′(xj) exp D ()
j=
One convenient definion of the sample points in x is
xj = (j − (N/))Δ, ()
where N is the number of sample points in the aperture. Note that this definion of the x-coordinate is equivalent to applying a coordinate transform equivalent to the “shi” operaon described in the lectures, and results in Fourier phases which are closer to zero compared to a more simple linear relaonship between xj and j.
. Tasks
Core Task : Write a program that will calculate the diffracon paern of a gen- eral -dimensional complex aperture in the far field of the aperture using FFT tech- niques. The program should calculate the intensity in the paern across the screen, which you should plot using the correct y coordinates (in metres or microns for ex- ample). Label your coordinates.
Test this program for the specific case of a slit in the centre of an otherwise blocked aperture: take the single slit to have width d in the centre of an aperture of total
extent L. For definiteness, use λ = nm, d microns, D = .m and L =mm. Overlay on your plot the theorecal value of the intensity paern expected.
Core Task : Now calculate and plot the Fraunhofer diffracon paern of a sinu- soidal phase grang. This grang is a slit of extent d = mm, outside of which the transmission is zero. Within |x| < d/, the transmission amplitude is ., and the phase of A is
φ(x) = (m/) sin(πx/s)
where s is the spacing of the phase maxima, and can be taken as microns for this problem. For this calculaon, use m = . The Fresnel distance d/λ is m, so calculate the paern on a screen at D = m. What do you noce about the resulng paern?
Supplementary Task : Now modify your program so that the calculaon is accu- rate even in the near-field by adding a phase correcon to the aperture funcon as defined by Equaon . Repeat your calculaons in the previous two tasks for D = mm for the slit, and D = .m for the phase grang, and plot the results. Do the intensity paerns look sensible?
. Hints
Recall the FFT definion:
Hj =
which maps N me-domain samples hm into N frequencies, which are
• j = is zero frequency.
• For ≤ n ≤ (N/), we have posive frequencies (j/N) × (/Δ).
• For (N/) + ≤ j ≤ (N − ) we have negave frequencies which we compute
as ((j/N) − ) × (/Δ). (Remember the sequences are periodic).
Of course in this case we don’t have me-domain samples, but we can sll use the
FFT to carry out the transform.
The complex aperture funcon will be represented by an array of N discrete com- plex values along the aperture, encoding the real and imaginary parts of A(x). Each complex value represents the aperture’s transmiance over a small length Δ of the aperture, so that NΔ is the total extent of the aperture.
N−
hmeπimj/N () fj= j ()
m=
NΔ
You can think of frequencies (j/N) × (/Δ), running from j = to (N − ), with
Choose appropriate values for N and Δ to make sure you can represent the whole aperture of maximum extent L well enough. Bear in mind that Fast Fourier Trans- form calculaons are fastest when the transform length is a power of , and that you want Δ to be small enough to resolve the features of the aperture. (Computers are fast these days though; so you can use small values of Δ and correspondingly large values of N. In pracce, for such a small problem, the use of N as a power of is not necessary, but it is important if performance is crical.)
For the slit problem you can use the funcon to set up an array of the appropriate size filled with zeros and then set the locaons where the slit is transparent by assigning non-zero values to “slices” e.g. .
You now need a roune to calculate the fast Fourier transform (FFT). You can use the funcon, which is straighorward to use (see the code ex- amples from the lectures). It accepts a real or a complex input and produces a complex output. You can call the FFT roune, write out the result to a text file as a and then read it into as separate program which uses for plong. You can do plong in the same program if you prefer. Compung and plong in the same program is not recommended if the computaon me is long (more than a few seconds), but this should not be the case here.
You will need to compute the intensity of the paern. You might also want to plot the aperture funcon to make sure you have calculated it correctly.
Now think carefully about the coordinates associated with your calculated paern. The discussion at the beginning of this secon reminds you about how the fre- quencies appear in the FFT’ed data. Can you understand the form of the intensity paern you have derived?
To plot the intensity on the screen as a funcon of actual distance y, you need to work out how to convert the pixels in the Fourier transform into distances on the screen y. To do this you first need to compare carefully Equaon and Equaon (also referring to Equaon ) which should tell you how to derive y at each pixel value. In addion, by interpreng the second half of the transform as negave frequencies (or y values in this case) you should be able to plot the intensity paern as a funcon of y for posive and negave y, and plot over this the matching sinc funcon for a slit. If you are having difficules with this step it may help to revise your notes from the lectures concerning the locaon of negave frequencies in the output of an FFT.