When collimated beam of light is passed through a slit, the beam diffracts and the diffraction pattern can be seen on a screen.

The resulting diffraction pattern can be derived to be of the form

$$ u(x,y,z)=\mathcal{F}^{-1}\left[A(k_x,k_y)e^{-iz\sqrt{k^2-k_x^2-k_y^2}}\right] $$

where $u(x,y,z)$ is the light distribution and $A(k_x, k_y)$ is the spatial frequency domain representation of $u$, and is given by

$$ A(k_x,k_y)=\mathcal{F}[u(x,y,0)]. $$

Note that $\mathcal{F}$ and $\mathcal{F^{-1}}$ is the Fourier and inverse Fourier transform respectively and $z=0$ is defined to be at the slit.

If you want to see the full derivation of this, see Don’t drink and derive → Optics → Fraunhofer Diffraction

Simulation Results

Python Code

import numpy as np
import scipy as sp
from scipy.fft import *
import imageio
import cv2

import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap
from matplotlib import animation 
from matplotlib.animation import PillowWriter
import pint # Deals with units

u = pint.UnitRegistry()

def foptics(shape='circle', a=0.1, d=0.15, bound=2):
    d = d * u.mm # distance between slits (only relevant for double slit)
    a = a * u.mm
    wavelength = 660 * u.nm

    x = np.linspace(-bound, bound, 2000) * u.mm
    X, Y = np.meshgrid(x, x)
		
    # Put slit shapes here!
    if shape == "doubleslit":
        u_0 = (np.abs(X-d/2) < a/2) * (np.abs(Y) < 0.5 * u.mm) + (np.abs(X + d/2) < a/2) * (np.abs(Y) < 0.5 * u.mm) # DOUBLE SLIT
    
    if shape == "singleslit":
        u_0 = (np.abs(X) < a/2) * (np.abs(Y) < 0.5 * u.mm) # SINGLE SLIT
    
    if shape == "circle":
        u_0 = X**2 + Y**2 < (a*2)**2 # CIRCULAR APERTURE
    
    if shape == "hexagon":
        hx = np.abs(X) < a
        hy = np.abs(Y) < np.sqrt(3) * a / 2
        hxy = np.abs(np.sqrt(3) * X - Y) < np.sqrt(3) * a
        hxy2 = np.abs(np.sqrt(3) * X + Y) < np.sqrt(3) * a
        u_0 = hx * hy * hxy * hxy2 # HEXAGONAL APERTURE
    
    if shape == "triangle":
        u_0 = (Y > -np.sqrt(3) * (X - a)) & (Y > np.sqrt(3) * (X + a)) & (Y < np.sqrt(3) * a) # EQUILATERAL TRIANGLE
    
    if shape == "squarearray":
        spacing = 2 * a
        u_0 = np.logical_or.reduce([(np.abs(X - i * spacing) < a/2) * (np.abs(Y - j * spacing) < a/2) for i in range(-2, 3) for j in range(-2, 3)]) # SQUARE ARRAY
    
    if shape == "ring":
        u_0 = ((X**2 + Y**2) < (2*a)**2) & ((X**2 + Y**2) > (a)**2) # RING SHAPE

    u_0 = u_0.astype(float)

    A = fft2(u_0)
    kx = fftfreq(len(x), np.diff(x)[0]) * 2 * np.pi
    Kx, Ky = np.meshgrid(kx, kx)

    def get_U(z, k):
        return ifft2(A*np.exp(1j * z * np.sqrt(k**2 - Kx**2 - Ky**2)))
    
    k = 2*np.pi/wavelength
    d = 3 * u.cm
    U = get_U(d, k)
    fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(24, 8))

    axes[0].pcolormesh(X, Y, u_0)
    axes[0].set_xlabel('x [mm]')
    axes[0].set_ylabel('y [mm]')
    axes[0].set_title('Slit Shape')

    axes[1].pcolormesh(X, Y, np.abs(U), cmap='inferno')
    axes[1].set_xlabel('$x$ [mm]')
    axes[1].set_ylabel('$y$ [mm]')
    axes[1].set_title('Diffraction Pattern on the Screen')

    axes[2].plot(x, np.abs(U)[250])
    axes[2].set_xlabel('$x$ [mm]')
    axes[2].set_ylabel('$u(x,y,z)$ [sqrt of intensity] [mm]')
    axes[2].set_title('Light Amplitude as a function of $x$')