\documentclass[11pt]{article}
\usepackage{amsmath}
\usepackage{autobreak}
\usepackage{lineno,hyperref}
\usepackage[table,x11names,dvipsnames,table]{xcolor}
\usepackage{authblk}
\usepackage{subcaption,booktabs}
\usepackage{graphicx}
\usepackage{multirow}
\usepackage{listings}
\usepackage{color}
\definecolor{blue}{rgb}{0.02,0.65,0.90}
\definecolor{dkgreen}{rgb}{0.25,0.63,0.17}
\definecolor{gray}{rgb}{0.5,0.5,0.5}
\definecolor{mauve}{rgb}{0.53,0.22,0.94}
\lstset{frame=tb,
  language=Python,
  aboveskip=3mm,
  belowskip=3mm,
  showstringspaces=false,
  columns=flexible,
  basicstyle={\small\ttfamily},
  numbers=none,
  numberstyle=\tiny\color{gray},
  keywordstyle=\color{blue},
  commentstyle=\color{dkgreen},
  stringstyle=\color{mauve},
  breaklines=true,
  breakatwhitespace=true,
  tabsize=3
}
\usepackage[nolist,nohyperlinks]{acronym}
\usepackage[superscript]{cite}
\usepackage{tabularx}
\usepackage{float}
\usepackage[group-separator={,}]{siunitx}
\usepackage{geometry}
 \geometry{
 a4paper,
 papersize={210mm,279mm},
 left=12.73mm,
 top=20.3mm,
 marginpar=3.53mm,
 textheight=238.4mm,
 right=12.73mm,
 }

\setlength{\columnsep}{6.54mm}

%\linenumbers %%% Turn on line numbers here

\renewcommand{\familydefault}{\sfdefault}

\captionsetup[figure]{labelfont=bf,textfont=normalfont}
%\captionsetup[subfigure]{labelfont=bf,textfont=normalfont}

%%%% comment out the below for the other title option
\makeatletter
\def\@maketitle{
\raggedright\newpage
  \noindent
  \vspace{0cm}
  \let\footnote\thanks{\hskip -0.4em \huge \textbf{{\@title}} \par}
    \vskip 1.5em
    {\large
      \lineskip.5em
      \begin{tabular}[t]{l}
      \raggedright\@author\end{tabular}\par}
    \vskip 1em
    \@date\par
  \vskip 1.5em
  }
\makeatother

\begin{document}

\title{Exercise 2 Report}

\author[1]{Ceres Milner}

\affil[1]{Department of Physics, University of Bristol}
\renewcommand\Affilfont{\itshape\small}% chktex 6 

\date{01.12.2025}
\maketitle
Word count: 1741
\begin{abstract}

	This exercise aimed to use numerical integration methods to investigate diffraction of light through an aperture onto a screen. Two calculation methods were used, a numerical calculation and Monte Carlo integration. The numerical method gave cleaner, less `noisy' data, but the Monte Carlo method was considerably quicker.  

\end{abstract}

\section{Introduction}

In this exersise we will use numerical methods to investigate the diffraction of light when passed through an aperture of different shapes. We will use two methods, first pure numerical integration methods, and then the Monte Carlo method, and compare the accuraccy and speed of these two methods. We will first create a 1 dimensional plot of the diffraction to ensure the results are as expected, and then 2 dimensional plots, using both square and curcular apertures.

\section{Theory and Methods}

To calculate the diffraction pattern of light passing through a single aperture, we will use the Fresnel diffraction integral, given by:
\begin{equation}
	E(x, y, z)=\frac{i^{ikz}}{i\lambda z}\int_{-\infty}^\infty\int_{-\infty}^\infty E(x^\prime,y^\prime)\text{exp}\bigg\{\frac{ik}{2z}\big[(x-x^\prime)^2+(y-y^\prime)\big]\bigg\}\mathrm{d}x^\prime\mathrm{d}y^\prime
\end{equation}
To allow for this to be calculated computtionally, we will simplify it, integrating over the area of the aperture, and seperating out the real and imaginary parts of the equation:
\begin{align}
	E(x,y,z)=\frac{kE_0}{2\pi z}\int_{x_1^\prime}^{x_2^\prime}\int_{y_1^\prime(x^\prime)}^{y_2^\prime(x^\prime)}\cos\bigg\{\frac{k}{2z}\Big[(x-x^\prime)^2+(y-y^\prime)^2\Big]\bigg\}\mathrm{d}x^\prime\mathrm{d}y^\prime \nonumber \\
	+i\frac{kE_0}{2\pi z}\int_{x_1^\prime}^{x_2^\prime}\int_{y_1^\prime(x^\prime)}^{y_2^\prime(x^\prime)}\sin\bigg\{\frac{k}{2z}\Big[(x-x^\prime)^2+(y-y^\prime)^2\Big]\bigg\}\mathrm{d}x^\prime\mathrm{d}y^\prime\label{eq:fresnel}
\end{align}

\section{Explanation of Code}

Initially, the functions for the kernels of the real and imaginary parts of equation~\ref{eq:fresnel} are set up, as well as setting some universal constants:
\begin{lstlisting}
def Fresnel2dreal(yp, xp, y, x, k, z): #Define functions for integral kernels
    kernel = np.cos((k/(2*z))*((x-xp)**2+(y-yp)**2))
    return kernel

def Fresnel2dimag(yp, xp, y, x, k, z):
    kernel = np.sin((k/(2*z))*((x-xp)**2+(y-yp)**2))
    return kernel

c = 3e8 #Universal constants
e0 = 8.85e-12
\end{lstlisting}
A function for part 1 is then defined:
\begin{lstlisting}
def plot1D(aperture, z, k, screen_range, resolution): #Function for part 1
    def genData(aperture, z, k, screen_range, resolution): #Function to generate data
        y = 0 #As in 1D

        xp1=yp1=-aperture/2 #Set range over which we are integrating
        xp2=yp2=aperture/2
    
        xs = np.linspace(-screen_range/2, screen_range/2, num=resolution) #Generate values to calculate for
        intensities = []

        constant = k/2*np.pi*z #Relative intensity constant
        
        for x in tqdm(xs): #Loop through all values of x
            realpart, realerror = integrate.dblquad(Fresnel2dreal, xp1, xp2, yp1, yp2, args=(y, x, k, z), epsabs=1e-10, epsrel=1e-10) #Calculate both real and imaginary parts
            imagpart, imagerror = integrate.dblquad(Fresnel2dimag, xp1, xp2, yp1, yp2, args=(y, x, k, z), epsabs=1e-10, epsrel=1e-10)

            I = c*e0*((realpart*constant)**2+(imagpart*constant)**2) #Combine parts and constants to get intensity
            intensities.append(I) #Add calculated intensity to list
        return xs, intensities

    ax = plt.axes()
    xs, intensities = genData(aperture, z, k, screen_range, resolution) #Plot intensity against distance
    ax.plot(xs, intensities)
    plt.xlabel("Position (m)")
    plt.ylabel("Relative Intensity")
    plt.title("1D diffraction")
    plt.show()
\end{lstlisting}
A sub function is defined with inputs of aperture, z, k screen range and resolution, to generate the data to plot. It begins by setting y=0 as it generates a 1D plot, and then setting the limits of the integration, defined by the aperture. An array of x values is then created, and these are the positions on the screen which we will calculate the integral for. The function then iterates through this array, and for each value of x it uses the dblquad function to calculate the real and imaginary part of the integral for the given x value. It then takes these results and combines them with the universal constants to give the intensity at that point, and then appends this value to a list of intensities. This genData function is then called and the data is saved to a variable, and the data is then plotted, creating a graph of intensity against position.

The function for part 2 is then defined:
\begin{lstlisting}
def plot2Drectangular(aperture, z, k, screen_range, resolution): #Function for part 2
    def genData(aperture, z, k, screen_range, resolution): #Function to generate data
        xp1=yp1=-aperture/2 #Set integration limits
        xp2=yp2=aperture/2
        
        xs = np.linspace(-screen_range/2, screen_range/2, num=resolution) #Generate values to integrate for
        ys = np.linspace(-screen_range/2, screen_range/2, num=resolution)

        intensities = []

        constant = k/2*np.pi*z #Relative intensity constant
        completion = 0

        for y in tqdm(ys):
            xIntensities = []
            for x in xs:
                realpart, realerror = integrate.dblquad(Fresnel2dreal, xp1, xp2, yp1, yp2, args=(y, x, k, z), epsabs=1e-10, epsrel=1e-10)#Calculate both parts
                imagpart, imagerror = integrate.dblquad(Fresnel2dimag, xp1, xp2, yp1, yp2, args=(y, x, k, z), epsabs=1e-10, epsrel=1e-10)

                I = c*e0*((realpart*constant)**2+(imagpart*constant)**2)#Combine both parts and constants
                xIntensities.append(I)
            intensities.append(xIntensities)
        intensities = np.array(intensities)
        return intensities

    intensity = genData(aperture, z, k, screen_range, resolution) #Generate the data
    extents = (-screen_range/2,screen_range/2,-screen_range/2,screen_range/2) #Set limits of the plot

    plt.imshow(intensity,vmin=0.0,vmax=1.0*intensity.max(),extent=extents,\
	origin="lower",cmap="nipy_spectral_r") #Plot the graphs
    plt.colorbar()
    plt.xlabel("X Position (m)")
    plt.ylabel("Y Position (m)")
    plt.title("2D diffraction with a rectangular aperture")
    plt.show()
\end{lstlisting}
The function for the second part is quite similar to the first part. Again, a function to generate the data is defined. It sets the linits for the integration, and this time generates an array of values for both x and y. It then loops through both the y and x lists, with the x loop nested inside the y loop. Each iteration it calculates the integral for the current location, and adds it to a 2D array representing the positions on the screen. Each row of pixels is a list, and the 2D array is a list of these rows. The genData function is then called and its result is saved to intensities, and the extents of the 2D plot are set to the size of the screen. The plot is then generated, with the scale set to be from 0 to the maximum value in the intensities array.

The function for part 3 is then defined:
\begin{lstlisting}
def plot2Dcircular(aperture, z, k, screen_range, resolution): #Function for part 3
    def genData(aperture, z, k, screen_range, resolution):
        xp1=-aperture/2 #Set integration limits
        xp2=aperture/2

        def yp1func(xp):
            return -np.sqrt((aperture/2)**2-(xp**2)) #Define y limits in terms of x

        def yp2func(xp):
            return np.sqrt((aperture/2)**2-(xp**2))

        xs = np.linspace(-screen_range/2, screen_range/2, num=resolution) #Generate values to integrate for
        ys = np.linspace(-screen_range/2, screen_range/2, num=resolution)

        intensities = []

        constant = k/2*np.pi*z

        for y in tqdm(ys):
            xIntensities = []
            for x in xs:
                realpart, realerror = integrate.dblquad(Fresnel2dreal, xp1, xp2, yp1func, yp2func, args=(y, x, k, z), epsabs=1e-10, epsrel=1e-10) #Calculate functions using circular limits
                imagpart, imagerror = integrate.dblquad(Fresnel2dimag, xp1, xp2, yp1func, yp2func, args=(y, x, k, z), epsabs=1e-10, epsrel=1e-10)

                I = c*e0*((realpart*constant)**2+(imagpart*constant)**2)
                xIntensities.append(I)
            intensities.append(xIntensities)
        intensities = np.array(intensities)
        return intensities

    intensity = genData(aperture, z, k, screen_range, resolution)
    extents = (-screen_range/2,screen_range/2,-screen_range/2,screen_range/2)
    
    plt.imshow(intensity,vmin=0.0,vmax=1.0*intensity.max(),extent=extents,\
	origin="lower",cmap="nipy_spectral_r")
    plt.colorbar()
    plt.xlabel("X Position (m)")
    plt.ylabel("Y Position (m)")
    plt.title("2D diffraction with a circular aperture")
    plt.show()
\end{lstlisting}
This operates very similarly to the part 2 function. The only difference is the addition of the yp1function and yp2function. These define the y prime limits of the integral in terms of the x prime value, and these are passed to the dblquad functions in order to calculate the integrals using a circular aperture.
The function for part 4 is then defined:
\begin{lstlisting}
def monte(aperture, z, k, screen_range, resolution, samples): #Function for part 4
    N = samples #Number of samples

    def doubleInteg(x, y, xp, yp, z, k, aperture): #Function to return calculated values for each sample
        values = []
        for i in range(len(xp)):
            if (xp[i]**2+yp[i]**2) > (aperture/2)**2: #Check if point is in the aperture
                values.append(0)
            else:
                value = np.exp(((1j*k)/(2*z))*((x-xp[i])**2+(y-yp[i])**2)) #Calculate value if in aperture
                values.append(value.imag)
        return np.array(values)
    
    def monteCarlo(x, y, z, k, aperture): #Perform onte carlo method
       xp = np.random.uniform(low=(-aperture/2), high=aperture/2, size=N) #Generate random samples 
       yp = np.random.uniform(low=(-aperture/2), high=aperture/2, size=N)
       values = doubleInteg(x, y, xp, yp, z, k , aperture) #Calculate value for each pair of samples
       mean = values.sum()/N
       meansq = (values*values).sum()/N
       integral = aperture*mean #Calculate the integral end error
       error = aperture*np.sqrt((meansq-mean*mean)/N)
       return integral, error
         
    def genData(aperture, z, k, resolution, screen_range): #Function to generate the data

        xs = np.linspace(-screen_range/2, screen_range/2, num=resolution) #Generate values to integrate for
        ys = np.linspace(-screen_range/2, screen_range/2, num=resolution)


        intensities = []
        
        constant = k/(2*np.pi*z)

        for y in tqdm(ys):
            xIntensities = []
            for x in tqdm(xs):
                integral, error = monteCarlo(x, y, z, k, aperture) #Find integral vaule using monte carlo method
                I = c*e0*constant*integral
                xIntensities.append(I)
            intensities.append(xIntensities)
        intensities = np.array(intensities)
        return intensities
    
    intensity = genData(aperture, z, k, resolution, screen_range) #Generate data
    extents = (-screen_range/2,screen_range/2,-screen_range/2,screen_range/2)

    for y in range(len(intensity)):
        for x in range(len(intensity[y])):
            if intensity[y][x] < 0.05*intensity.max():
                intensity[y][x] = 0

    plt.imshow(intensity,vmin=0,vmax=1.0*intensity.max(),extent=extents,\
	origin="lower",cmap="nipy_spectral_r")
    plt.colorbar()
    plt.xlabel("X Position (m)")
    plt.ylabel("Y Position (m)")
    plt.title("2D diffraction through Monte Carlo")
    plt.show()
\end{lstlisting}
The function doubleInteg is defined. This takes values of x, y, z, k and aperture, as well as a list of random samples of xp and yp. For each value pair of xp and yp, it test if the coordinate is in the circular aperture. If it is, the kernel function is calculated and appended to the array of values, and of the coordinate is outside of the aperture, the value is said to be 0 and this is appended to the list of values. The monteCarlo function is then defined. This first generates N random samples of xp and yp in two arrays. It then calls the doubleInteg function, passing these arrays of samples as well as the other needed values, and saves the resulting list of values to the values variable. From this list of values, it calculates both the mean and mean squared, and uses these to calculate an approximation of the integral, as well as the error associated with this integral, and returns both. As with the previous sections, the genData function is the defined, first creating arrays of x and y values, then looping through these arrays. This time however, instead of using the dblquad function it uses the monteCarlo function to calculate the integral, then calculates the intensity from this. Again, it saves each value to a 2D array, and returns this array. The genData function is then called and the array is saved to intensity, and the plot extents are set. As the Monte Carlo method generates a lot of low-level noise, the intensity is then looped through and all values less than 5\% of the maximum are set to 0. While this may remove some actual data, it greatly increases the signal to noise ration, allowing for the diffraction pattern to be observed more easily. The data is then plotted in the same manner as parts 2 and 3.

\section{Results and Discussion}

Throughout these results a wavelength of 1e-6m will be used. For far field diffraction, an aperture of 2e-5m and a z of 0.05m will be used, and for near field diffraction an aperture of 2e-4 and a z of 0.005m will be used.
\subsection{Section 1: 1D diffraction}
The 1D diffraction performs very well for far field diffraction, giving a graph in the expected shape, as shown in figure~\ref{fig:1d-far-field}
\begin{figure}
    \includegraphics[width=1\linewidth]{Images/1d-far-field.png}
    \caption{1D Far-field diffraction}\label{fig:1d-far-field}
\end{figure}
For near field diffraction, initially the plot was very noisy, with many spikes as shown in figure~\ref{fig:1d-far-field-wrong}. 
\begin{figure}
    \includegraphics[width=1\linewidth]{Images/1d-near-field-wrong.png}
    \caption{1D Near-field diffraction with incorrect tolerences}\label{fig:1d-far-field-wrong}
\end{figure}
To correct for this, the epsabs and epsrel arguments to 1e-10, which produced the correct diffraction pattern, as shown in figure~\ref{fig:1d-near-field}
\begin{figure}
    \includegraphics[width=1\linewidth]{Images/1d-near-field.png}
    \caption{1D Near-field diffraction}\label{fig:1d-near-field}
\end{figure}
\subsection{Section 2: 2D diffraction with a square aperture}
Section 2 of the code also behaves as expected. It begins to become quite time intensive at high resolutions due to the large number of integrations needed to be calculated. It produces the expected pattern for both near and far feild diffraction, as shown in figures~\ref{fig:2d-far-field} and~\ref{fig:2d-near-field}
\begin{figure}
    \includegraphics[width=1\linewidth]{Images/2d-far-field.png}
    \caption{2D Far-field diffraction}\label{fig:2d-far-field}
\end{figure}
\begin{figure}
    \includegraphics[width=1\linewidth]{Images/2d-near-field.png}
    \caption{2D Near-field diffraction}\label{fig:2d-near-field}
\end{figure}
As the near field diffraction takes longer to compute than the far field, the near field plot was generated using a lower resolution of 100$\times$100, whereas the far field plot was 400$\times$400.
\subsection{Section 3: 2D diffraction with a circular aperture}
Again, section 3 gave the expected results, as shown in figures~\ref{fig:2d-circular-far-field} and~\ref{fig:2d-circular-near-field}, however for this section the near field diffraction took much longer to compute than the far field, so again it was performed at a lower resolution. 
\begin{figure}
    \includegraphics[width=1\linewidth]{Images/2d-circular-far-field.png}
    \caption{2D Far-field diffraction with a circular aperture}\label{fig:2d-circular-far-field}
\end{figure}
\begin{figure}
    \includegraphics[width=1\linewidth]{Images/2d-circular-near-field.png}
    \caption{2D Far-field diffraction with a circular aperture}\label{fig:2d-circular-near-field}
\end{figure}
\subsection{Section 4: 2D diffraction through Monte Carlo}
The Monte Carlo method computed much faster than the numerical integration. At low resolutions and a low number of samples, a weak pattern was observed in the centre of the plot, but the rest of the data was effectively just noise, as shown in figure~\ref{fig:monte-far-field-low-res-low-samples}. While this is very quick to compute, it is very hard to draw any data from it.
\begin{figure}
    \includegraphics[width=1\linewidth]{Images/monte-far-field-low-res-low-samples.png}
    \caption{Monte Carlo plot with low resolution and a low number of samples}\label{fig:monte-far-field-low-res-low-samples}
\end{figure}
When given a low resolution of 50$\times$50 pixels, but a high sample rate of 5000 samples, the plot was still preduced quite quickly, but now showed a strong pattern, shown in figure~\ref{fig:monte-far-field-low-res-high-samples}, although with this low resolution it is still hard to confirm the accuraccy of the pattern.
\begin{figure}
    \includegraphics[width=1\linewidth]{Images/monte-far-field-low-res-high-samples.png}
    \caption{Monte Carlo plot with low resolution and a high number of samples}\label{fig:monte-far-field-low-res-high-samples}
\end{figure}
For a plot with a high resolution of 500$\times$500 and 10 samples per pixel, the pattern was more distinguishable than the low resolution plot, but outside the middle of the plot the data quickly became noise and it is hard to make out any pattern, as shown in figure~\ref{fig:monte-far-field-high-res-low-samples}
\begin{figure}
    \includegraphics[width=1\linewidth]{Images/monte-far-field-high-res-low-samples.png}
    \caption{Monte Carlo plot with high resolution and a low number of samples}\label{fig:monte-far-field-high-res-low-samples}
\end{figure}
For a plot with both a high resolution of 600$\times$600 pixels and a high number of samples, 1000, a clear pattern is created, show in fgure~\ref{fig:monte-far-field-high-res-high-samples}, although this takes much longer to compute than the previous plots. This plot also contains regular circular patterns not present in the numerical integrations. When compared with the numerical method, the MonteCarlo method shows the same pattern, although with more noise present and some irregularities.
\begin{figure}
    \includegraphics[width=1\linewidth]{Images/monte-far-field-high-res-high-samples.png}
    \caption{Monte Carlo plot with high resolution and a high number of samples}\label{fig:monte-far-field-high-res-high-samples}
\end{figure}
Unlike the previous sections, the near field pattern does not take much longer than the far field, likely as the dblquad function tests for convergence, which takes longer with the integrations present in the nearfield situation, however this is not the case for the Monte Carlo method, so the near field situation takes aproximately the same time as the far firld situation. When using a high resolution of 400$\times$400 pixels and a high number of scamples of 1000, the Monte Carlo method also produces a clear pattern for the near field situation, shown in figure~\ref{fig:monte-near-field-high-res-high-samples}, and again this matches the pattern produced by the numerical integration
\begin{figure}
    \includegraphics[width=1\linewidth]{Images/monte-near-field-high-res-high-samples.png}
    \caption{Monte Carlo plot of near field diffraction with high resolution and a high number of samples}\label{fig:monte-near-field-high-res-high-samples}
\end{figure}
\subsection{Section 5: Time for data generation}
To compare the efficiency of the numerical method and the Monte Carlo method, plots were generated for a number of different resolutions using the numerical method and the Monte Carlo method at a nu,ber of different sample rates. A plot of time against resolution was then created. As the number of integrations performed is proportional to resolution$^2$, we would expect the time taken to also be proportional to resolution$^2$. As such, time was plotted against resolution$^2$ to produce straight lines, shown in figure~\ref{fig:times}. As expected, each line follows a very strong linear relationship, showing that time is indeed proportional to resolution$^2$. Additionally, we can see that the samples needed for the Monte Carlo method to match the numerical method in efficiency is approximately 1600, however even at theis number of samples the Monte Carlo method still gives less clear data than the numerical method. As such, to balance if the accuracy of a Monte Carlo plot with samples $\leq$ 1600 is enough, the Monte Carlo method should be used as it is much faster, but if a greater accuracy is needed, the numerical method should be used.
\begin{figure}
    \includegraphics[width=1\linewidth]{Images/times.png}
    \caption{A plot of time against resolution$^2$ for several different plotting methods}\label{fig:times}
\end{figure}

\section{Conclusion}

All numerical methods gave the expected results, hadeling both near and far field diffraction, however at higher resolutions the 2 dimensional plots began to take a large amount of time to compute. Using the Monte carlo method, if either the resolution of number of samples was too low, the produced plot gave little useful data, however with sufficiently large resolution and number of samples, it gave a good plot whilst still being more efficient than the numerical method.

\end{document}
