PP/JAI Web>Simulation>PythonOptics>PythonOpticsGemmaSmithDiffraction (21 Nov 2013, RobAinsworth)

Fourier Optics - 'Diffraction calculations '

Aims

An ability to simulate any optical system
Compile a library of optical functions
Gain an understanding of Python
Learn about Frauhofer and Fresnel integrals

Background

There are some basic pieces of information that are need in this project. There are as follows:

Fourier Transform

This is the general form of a Fourier transform for 2-D which is needed so that in complex integrals found in Fresnel can be transferred into a more useful form

$\begin{displaymath} G(x,y) = \int\int_{-\infty}^\infty g(\xi,\eta) e^{-i2\pi(\xi x + \eta y)} d\xi \end{displaymath}$

(1)

Discrete Fourier Transform

This is the general form of a Fourier transform for 2-D which is the equation a computer uses to do a Fourier transform. A computer cant do a calculation over an infinite space, it and only deal with single points at a time and then summing them, so as can be seen in the equation in stead of dealing with x and y, the computer deals with integers p and q.

$\begin{displaymath} DFT(g_{m,n}) = G_{p,q} = \sum_{m=0}^{M-1}\sum_{n=0}^{N-1} g_{n,n}e^{-2\pi(\frac{mp}{M} + \frac{nq}{N})} \end{displaymath}$

(2)

Functions

There are three basic function that are useful:

Rectangular

This is used to simulate a rectangular aperture

$\begin{displaymath} \begin{equation*} rect(x) = \left\{ \begin{array}{rl} 1 & |x| < \frac{1}{2}\\ \frac{1}{2} & |x| = \frac{1}{2}\\ 0 & \text{otherwise } \end{array} \right. \end{equation*} \end{displaymath}$

(3)

2-D Plot of Rectangular Function:

rectangular.py.txt: Rectangluar Function

Circular

This is used to simulate a circular aperture

$\begin{displaymath} \begin{equation*} circ(r) = \left\{ \begin{array}{rl} 1 & r < 1\\ \frac{1}{2} & r = 1\\ 0 & \text{otherwise } \end{array} \right. \end{equation*} \end{displaymath}$

(4)

where

$\begin{displaymath} r = \sqrt{x^2 + y2} \end{displaymath}$

(5)

2-D Plot of Circular Function:

circular.py.txt: Circular Function

Gaussian

This is the most useful for real case simulations as most lasers would have a Gaussian wave form.

$\begin{displaymath} G(x,y) = e^{- \frac{x^2 + y^2}{2 w^2}} \end{displaymath}$

(6)

2-D Plot of Gaussian Function:

Gaussian.py.txt: Gaussian Function

Book work

Huygens-Fresnel Principle

All this project stems from the Huygen's Fresnel Principle which describes what happens to a wavefield when it propagates some perpendicular distance in z between to parallel planes. The equation for this is as follows

$\begin{displaymath} U(x,y) = \frac{z}{i\lambda} \int\int_{\Sigma} U(P_1)\frac{e^{ikr_{01}}}{r_{01}^2} ds \end{displaymath}$

(7)

where

$\begin{displaymath} r_{01} = \sqrt{z^2 + (x - \xi )^2 + ( y - \eta )^2} \end{displaymath}$

(8)

From this two approximation are applied so as to enable a relationship to a fourier transform to be made.

The Fresnel Approximation

With Fresnel a binomial expansion is applied to r.

The general form of the expansion is

$\begin{displaymath}\ \sqrt{1+b} = 1 + \frac{1}{2} b + \frac{1}{8} b^2 + ... \end{displaymath}$

(9)

Now we rearrange r by removing z fromt he square root as to match up with the general expresion

$\begin{displaymath} r_{01} = z \sqrt{1 + (\frac{x - \xi}{z})^2 + (\frac{y - \eta}{z})^2} \end{displaymath}$

(10)

So now we can presides with the expansion, taking only the first 3 terms as relevant

$\begin{displaymath} r_{01} \approx z [1 + \frac{1}{2}(\frac{x - \xi}{z})^2 + \frac{1}{2}(\frac{y - \eta}{z})^2] \end{displaymath}$

(11)

Now this new r is substituted into huygens-Fresnel, qith the exponinent this is a simple thing but int he fraction r the contribution from the terms involving x and y are so small that they become in consequential. all this leads to the following equation which is the general form of the Fresnel approximation

$\begin{displaymath} U(x,y) = \frac{e^{ikz}}{i\lambda z} \int\int_{-\infty}^\infty U(\xi,\eta) exp(\frac{ik[(x-\xi)^2+(y-\eta)^2]}{2z}) d\xi d\eta \end{displaymath}$

(12)

The Fraunhofer Approximation

With Fraunhofer the Fresnel approximation is taken and generalised for very large distances where

$\begin{displaymath} z >> \frac{k(\xi^2 + \eta^2)_{max}}{2} \end{displaymath}$

(13)

This means that this part of the expanded exponent can be ignored to give

$\begin{displaymath} U(x,y) = \frac{e^{ikz}e^\frac{ik(x^2 + y^2)}{2z}}{i\lambda z} \int\int_{-\infty}^\infty U(\xi, \eta)exp(\frac{-ik(\xi x + \eta y)}{z}) d\xi d\eta \end{displaymath}$

(14)

The condition which should be satisfied so that this approximation can be used is

$\begin{displaymath} z > \frac{2D^2}{\lambda} \end{displaymath}$

(15)

where D is the linear dimension of the aperture.

Book to Computer

How a human and computer approach a problem is very different. So all the analytical calculation that have been done up to this point must now be transformed into the discrete calculation a computer will preform.

Analytical Vs. Discrete

A Fourier transform is normally an analytical calculation, but a computer does it discretely gemma smells, this two approaches can not be exactly related, therefore a correction factor must be found to enable an exact comparison and further more enable the computer to do the right calculation. this is done as follows.

Comparing Fourier sections

$\begin{displaymath} \int\int_{-\infty}^\infty U(\xi,\eta)e^{\frac{-ik(\xi x + \eta y)}{z}} d\xi d\eta = \sum_{m=0}^{M-1}\sum_{n=0}^{N-1} g_{n,n}e^{-2i\pi(\frac{mp}{M} + \frac{nq}{N})} \end{displaymath}$

(16)

equationing releven parts

$\begin{displaymath} U(\xi,\eta) = g_{n,n} \end{displaymath}$

(17)

$\begin{displaymath} e^{\frac{-ik(\xi x + \eta y)}{z}} = e^{-2i\pi(\frac{mp}{M} + \frac{nq}{N})} \end{displaymath}$

(18)

log both sides

$\begin{displaymath} \frac{-ik(\xi x + \eta y)}{z} = -2i\pi(\frac{mp}{M} + \frac{nq}{N}) \end{displaymath}$

(19)

replacing k

$\begin{displaymath} {\frac{-i2\pi(\xi x + \eta y)}{\lambda z}} &=& -2i\pi(\frac{mp}{M} + \frac{nq}{N}) \end{displaymath}$

(20)

cancel down

$\begin{displaymath} {\frac{(\xi x + \eta y)}{\lambda z}} = \frac{mp}{M} + \frac{nq}{N} \end{displaymath}$

(21)

equating for:

- x which is related to p

$\begin{displaymath} \frac{\xi x}{\lambda z} = \frac{mp}{M} \end{displaymath}$

(22)

rearranged

$\begin{displaymath} x = \frac{mp\lambda z}{\xi M} \end{displaymath}$

(23)

- y which is related to p

$\begin{displaymath} \frac{\eta y}{\lambda z} = \frac{nq}{N} \end{displaymath}$

(24)

rearranged

$\begin{displaymath} y = \frac{nq\lambda z}{N\eta} \end{displaymath}$

(25)

we also know that

$\begin{displaymath} x = p d_x \end{displaymath}$

(26)

$\begin{displaymath} \xi = m d_{\xi} \end{displaymath}$

(27)

$\begin{displaymath} \frac{p d_x m d_{\xi}}{\lambda z} = \frac{mp}{M} \end{displaymath}$

(28)

canceling m and p

$\begin{displaymath} \frac{d_x d_{\xi}}{\lambda z} = \frac{1}{M} \end{displaymath}$

(29)

Rearanging for dx

$\begin{displaymath} d_x = \frac{\lambda z}{M d_{\xi}} \end{displaymath}$

(30)

Similarly for dy

$\begin{displaymath} d_y = \frac{\lambda z}{N d_{\eta}} \end{displaymath}$

(31)

General Programing

There are element that will be constantly reapering through out this equation, i.e.

Inputs
- wavelength - wl
- distance from aperture - z
- Size of Aperture - w
- Spacing between points - d

Output
- Fourier Transform - ft
- New Wave field - U
- Intensity - U*U.cong()

1-D Rectangular Apertures

To enable a easy understand when starting to program, we start with doing this calculation in 1-D, this also means that the calculation is easier to understand.

Fraunhofer

We start with the Founhofer approximation as its the simpler of the two.

This is the start wave field put in:

$\begin{displaymath} U(\xi) = rect(\frac{\xi}{2w}) \end{displaymath}$

(32)

This is the end intensity when the calculation is done by hand:

$\begin{displaymath} I(x) = \frac{A^2}{\lambda^2 z^2} sinc^2(\frac{2wx}{\lambda z}) \end{displaymath}$

(33)

where

$\begin{displaymath} A = 2w \end{displaymath}$

(34)

1-D Plot for a Franuhofer Calculation, with: *wavelength = 60mm *z = 4km *w = 5m:

fraunhoferRect.py.txt: Code for Fraunhofer Rectangular Apperture

Fresnel

Now the complex Fresnel approximation is done, this is the same as the Franhofer but with the addition of an extra phase within the Fourier transform.

This is the start wave field put in:

$\begin{displaymath} U(\xi) = rect(\frac{\xi}{2w}) \end{displaymath}$

(35)

This is the end intensity when the calculation is done by hand:

$\begin{displaymath} I(x) = \frac{1}{2 \sqrt{\lambda z}} \{ [(C(\alpha 2) - C(\alpha 1)]^2 + [S(\alpha 2) - S(\alpha 1)]^2\} \end{displaymath}$

(36)

where we use the Fresnel integrals

$\begin{displaymath} C(z) = \int_{0}^{z} \cos(\frac{\pi t^2}{2}) dt \end{displaymath}$

(37)

$\begin{displaymath} S(z) = \int_{0}^{z} \sin(\frac{\pi t^2}{2}) dt \end{displaymath}$

(38)

and where

$\begin{displaymath} \alpha 1 = - \sqrt{\frac{2}{\lambda z}} (w + x) \end{displaymath}$

(39)

$\begin{displaymath} \alpha 2 = \sqrt{\frac{2}{\lambda z}}(w - x) \end{displaymath}$

(40)

Another interesting piece of information is the Fresnel Number:

$\begin{displaymath} N_f = \frac{w^2}{\lambda z} \end{displaymath}$

(41)

As can be seen in the example below as the Fresnel number increases the distance decreases, also it can be seen how as the distance gets closer to the one used for a Fraunhofer approximation the shapes become the same.For this example all the wavelengths are 60mm and the widths is 5m.

1-D Plot for a Fresnel Calculation, with: *z = 4km *Fresnel Number = 0.1:

1-D Plot for a Fresnel Calculation, with: *z = 400m *Fresnel Number = 1:

1-D Plot for a Fresnel Calculation, with: *z = 100m *Fresnel Number = 4:

1-D Plot for a Fresnel Calculation, with: *z = 40m *Fresnel Number = 10:

* fresnelRect.py.txt: Code for Fresnel Rectangular Apperture

2-D

The transformation from 1-D to 2-D is simple for the equation as all that is needed is to add in the y componts along with the x, it become more complex for the decrete points being used as now instead of just a line of points, it is necessary to create a grid of points. This is done by taking 2 arrays, x and y, and meshing them into a grid.

Classes

Up until this point all the code has been written as a list of commands. By transfering this code into classes it means that there will be less errors, as once a piece of code it correct it can then just be called on. Also it means that later when multiple elements are simulated simultaneously, this can be done easily by calling the relevent piece of code, meaning hundreds of lines of code can be compressed into just a few lines.

In this case the code has been split into two classes

Intensity - this is where all the information about a wave field, its grids and all plots, this means that no information about the intensity is lost at any point through the calculation.

intensity.py.txt: Intensity Class Code

Progation - this is a function class where all the transformation on a wave field are done.

propagate.py.txt: Propagation Class Code

Examples

These are examples of 2-D plots using the new classes system

2-D Plot of Rectangular Function through Fraunhofer at z = 4m and w = 0.5mm:

2-D Plot of Rectangular Function through Fresnel at z = 0.01m and w = 0.5mm:

2-D Plot of Circular Function through Fraunhofer at z = 0.1m and w = 0.05mm:

2-D Plot of Circular Function through Fresnel at z = 0.1, and w = 0.5mm:

2-D Plot of Gaussian Function through Fraunhofer at z = 20m and w = 0.05mm:

2-D Plot of Gaussian Function through Fresnel at z = 0.02m and w = 0.5mm:

rum.py.txt: Test Code using Classes system and in 2-D

Lens

When adding a lens it is just a case of ajusting the phase of the wave field. This is done by multipying the wave feild by the a phase factor before doing a normal fresnel approxiamtion.

Lens phase alteration

$\begin{displaymath} U_L = U_I e^{\frac{-ik}{2f}(x^2 + y^2)} \end{displaymath}$

(42)

Where f is the focal length

Example

This is an example of a gaussian wave field propagating thought a lens and then over sum distant to the focal point. With a focal length of 0.5m and w = 1mm

Plots just after lens:

These are intensity plot along z. NB - this distance in z is at the top of each graph:

This is a plot of the second moment, or in other words its a plot of the outside edge of constructive interference:

Lens.py.txt: Code for Lens

NB The lens is still having problems and therefore still being worked on

Latex rendering error!! dvi file was not created.

Attachments

Topic attachments
I	Attachment	History	Action	Size	Date	Who	Comment
png	2M.png	r1	manage	28.7 K	11 Aug 2009 - 16:07	GemmaSmith	This is a plot of the second moment, or in other words its a plot of the outside edge of constuctive interference
png	Ciri.png	r1	manage	53.8 K	07 Aug 2009 - 16:02	GemmaSmith	2-D Plot of Circular Function
png	Gau.png	r1	manage	61.0 K	07 Aug 2009 - 16:03	GemmaSmith	2-D Plot of Gaussian Function
txt	Gaussian.py.txt	r1	manage	0.1 K	07 Aug 2009 - 16:07	GemmaSmith	Gaussian Function
txt	Lens.py.txt	r1	manage	1.1 K	11 Aug 2009 - 16:48	GemmaSmith	Code for Lens
png	Rect.png	r1	manage	46.9 K	07 Aug 2009 - 16:05	GemmaSmith	2-D Plot of Rectangular Function
txt	circular.py.txt	r1	manage	0.2 K	07 Aug 2009 - 15:50	GemmaSmith	Circular Function
png	fra.png	r1	manage	43.0 K	07 Aug 2009 - 16:45	GemmaSmith	1-D Plot for a Franuhofer Calculation, with: wavelength = 60mm z = 4km *w = 5m
txt	fraunhoferRect.py.txt	r1	manage	0.8 K	07 Aug 2009 - 17:02	GemmaSmith	Code for Fraunhofer Rectangular Apperture
png	fre0.1.png	r1	manage	55.2 K	07 Aug 2009 - 16:58	GemmaSmith	1-D Plot for a Fresnel Calculation, with: z = 4km Fresnel Number = 0.1
png	fre1.png	r1	manage	50.0 K	07 Aug 2009 - 16:48	GemmaSmith	1-D Plot for a Fresnel Calculation, with: z = 400m Fresnel Number = 1
png	fre10.png	r1	manage	72.1 K	07 Aug 2009 - 16:59	GemmaSmith	1-D Plot for a Fresnel Calculation, with: z = 40m Fresnel Number = 10
png	fre4.png	r1	manage	57.7 K	07 Aug 2009 - 17:00	GemmaSmith	1-D Plot for a Fresnel Calculation, with: z = 100m Fresnel Number = 4
txt	fresnelRect.py.txt	r1	manage	1.4 K	07 Aug 2009 - 17:03	GemmaSmith	Code for Fresnel Rectangular Apperture
png	gua_fra.PNG	r1	manage	15.4 K	10 Aug 2009 - 17:00	GemmaSmith	2-D Plot of gaussian through fraunhofer
png	image.png	r1	manage	207.4 K	11 Aug 2009 - 16:01	GemmaSmith	Plots just after lens
png	intensity.png	r1	manage	102.8 K	11 Aug 2009 - 16:04	GemmaSmith	This are intensity plot along z NB - this distance in z is at the top of each graph
txt	intensity.py.txt	r1	manage	3.7 K	11 Aug 2009 - 14:21	GemmaSmith	Intensity Class Code
txt	propagate.py.txt	r1	manage	1.3 K	11 Aug 2009 - 14:21	GemmaSmith	Propagation Class Code
txt	rectangular.py.txt	r1	manage	0.2 K	07 Aug 2009 - 15:51	GemmaSmith	Rectangluar Function
txt	rum.py.txt	r1	manage	1.3 K	11 Aug 2009 - 14:11	GemmaSmith	Test Code using Classes system and in 2-D
png	testcirfra.png	r1	manage	31.6 K	11 Aug 2009 - 13:55	GemmaSmith	2-D Plot of Circular Function through Fraunhofer
png	testcirfre.png	r1	manage	223.8 K	11 Aug 2009 - 13:52	GemmaSmith	2-D Plot of Circular Function through Fresnel
png	testfreRec.png	r1	manage	330.4 K	11 Aug 2009 - 15:00	GemmaSmith	2-D Plot of Rectangular Function through Fresnel
png	testgaufra.png	r1	manage	27.4 K	11 Aug 2009 - 14:02	GemmaSmith	2-D Plot of Gaussian Function through Fraunhofer
png	testgaufre.png	r1	manage	307.0 K	11 Aug 2009 - 14:03	GemmaSmith	2-D Plot of Gaussian Function through Fresnel
png	testrecfra.png	r1	manage	35.1 K	11 Aug 2009 - 13:56	GemmaSmith	2-D Plot of Rectangular Function through Fraunhofer

Physics Webpages • RHUL Webpages • Campus Connect • Royal Holloway, University of London, Egham, Surrey TW20 0EX; Tel/Fax +44 (0)1784 434455/437520

Topic revision: r14 - 21 Nov 2013 - RobAinsworth

Copyright © 2008-2025 by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
Ideas, requests, problems regarding RHUL Physics Department TWiki? Send feedback