latex/schrodinger_simulation.tex

@@ -82,7 +82,7 @@
 
 \begin{document}
 
-\title{Simulating the Schrödinger wave equation using the Crank-Nicolson method in 2+1 dimensions}  % self-explanatory
+\title{Simulating the Schrödinger wave equation using the Crank-Nicolson scheme in 2+1 dimensions}  % self-explanatory
 \author{Cory Alexander Balaton \& Janita Ovidie Sandtrøen Willumsen \\ \faGithub \, \url{https://github.uio.no/FYS3150-G2-2023/Project-4}} % self-explanatory
 \date{\today}                             % self-explanatory
 \noaffiliation                            % ignore this, but keep it.
@@ -122,12 +122,14 @@
 
 \end{document}
 
-% Methods
+% Abstract: OK
-% P1: Theory, imag i = i, index i, j = \hat{i}, \hat{j}
+% Introduction: 
-% P2:
+% Methods: Wave packet, Heisenberg, Hyugen
-% 
+% Results: 
+% Conclusion: OK
+% Appendices: 
 
 % Results
-% P7:
+% P7: OK
 % P8:
 % P9:

latex/sections/abstract.tex

@@ -3,16 +3,18 @@
 \begin{document}
 \begin{abstract}
 We have simulated the two-dimensional time-dependent Schrödinger equation, to study 
-variations of the double-slit experiment. To solve the partial differential equations 
+variations of the double-slit experiment. To derive a discretized equation 
-we have applied the Crank-Nicolson scheme in 2+1 dimensions, to derive a discretized 
+we applied the Crank-Nicolson scheme in 2+1 dimensions. In addition, we have used 
-equation. In addition, we have used Dirichlet boundary conditions to express the 
+Dirichlet boundary conditions to express the equation in matrix form, and solve 
-equation in matrix form and solve it using the sparse matrix solver \verb|superlu|. 
+it using the sparse matrix solver \verb|superlu|. Our implementation, and choice 
-Our implementation, and choice of solver method, resulted in conserved total 
+of solver method, resulted in a deviation from conserved total probability on the 
-probability $\sum_{\ivec, \jvec} p_{\ivec, \jvec}^{n}=1$ for both the single and 
+scale $10^{-14}$, for both the single and double slit setup. To illustrate the time 
-double slit setup. To illustrate the time evolution of the probability function, 
+evolution of the probability function, we created colormap plots for time steps 
-we created colormap plots at time steps $t = [0, 0.001, 0.002]$. We also included 
+$t = [0, 0.001, 0.002]$. We also included separate plots for each time step of 
-separate plots for each time step of Re$(u_{\ivec, \jvec})$ and Im$(u_{\ivec, \jvec})$.
+Re$(u_{\ivec, \jvec})$ and Im$(u_{\ivec, \jvec})$. In addition, we determined the 
-In addition, we determined the normalized particle detection probability $p(y \ | \ x=0.8, t=0.002)$, 
+normalized particle detection probability $p(y \ | \ x=0.8, t=0.002)$, for single-, 
-for single-, double- and triple-slit.
+double- and triple-slit.
 \end{abstract}
 \end{document}
+
+% $| \sum_{\ivec, \jvec} p_{\ivec, \jvec}^{n} - 1 | \approx $ 

latex/sections/appendices.tex

@@ -25,10 +25,10 @@ We need the first derivative in respect to both time and position, as well as th
 Schrödinger contain $i$ at the lhs, factor it as
 \begin{align*}
     \frac{u_{\ivec, \jvec}^{n+1} - u_{\ivec, \jvec}^{n}}{\Delta t} &= \frac{1}{2i} \bigg[ F_{\ivec, \jvec}^{n+1} + F_{\ivec, \jvec}^{n} \bigg] \\
-    &= -\frac{i}{2} \bigg[ F_{\ivec, \jvec}^{n+1} + F_{\ivec, \jvec}^{n} \bigg] & \text{, where $\frac{1}{i} = -i$}
+    &= -\frac{i}{2} \bigg[ F_{\ivec, \jvec}^{n+1} + F_{\ivec, \jvec}^{n} \bigg] \ ,
 \end{align*}
-
+where $1/i = -i$. Using Equation \eqref{eq:schrodinger_dimensionless}, which is 
-Using Equation \eqref{eq:schrodinger_dimensionless}, we get
+found in Section \ref{ssec:schrodinger}, we get
 \begin{align*}
     u_{\ivec, \jvec}^{n+1} - u_{\ivec, \jvec}^{n} & -\frac{i \Delta t}{2} \bigg[ F_{\ivec, \jvec}^{n+1} + F_{\ivec, \jvec}^{n} \bigg] \\
     &= -\frac{i \Delta t}{2} \bigg[ - \frac{u_{\ivec+1, \jvec}^{n+1} - 2u_{\ivec, \jvec}^{n+1} + u_{\ivec-1, \jvec}^{n+1}}{2 \Delta x^{2}} \\
@@ -94,4 +94,55 @@ $(M-2)^{2} \times (M-2)^{2} = 9 \times 9$ given by
     0   &  0   &  0   &  0    &   0   &  r    &  0   &  r    &  b_{8}  \\
     \end{bmatrix}
 \end{align*}
+
+
+\section{Figures}\label{ap:figures}
+\begin{figure*}
+    \centering
+    \begin{subfigure}[b]{0.3\textwidth}
+        \centering
+        \includegraphics[width=\textwidth]{images/color_map_0_real.pdf}
+        \caption{Re$(u_{\ivec, \jvec})$ at time $t=0$.}
+        \label{fig:colormap_0_real}
+    \end{subfigure}
+    \hfill
+    \begin{subfigure}[b]{0.3\textwidth}
+        \centering
+        \includegraphics[width=\textwidth]{images/color_map_1_real.pdf}
+        \caption{Re$(u_{\ivec, \jvec})$ at time $t=0.001$.}
+        \label{fig:colormap_1_real}
+    \end{subfigure}
+    \hfill
+    \begin{subfigure}[b]{0.3\textwidth}
+        \centering
+        \includegraphics[width=\textwidth]{images/color_map_2_real.pdf}
+        \caption{Re$(u_{\ivec, \jvec})$ at time $t=0.002$.}
+        \label{fig:colormap_2_real}
+    \end{subfigure}
+
+    \newline
+
+    \begin{subfigure}[b]{0.3\textwidth}
+        \centering
+        \includegraphics[width=\textwidth]{images/color_map_0_imag.pdf}
+        \caption{Im$(u_{\ivec, \jvec})$ at time $t=0$.}
+        \label{fig:colormap_0_imag}
+    \end{subfigure}
+    \hfill
+    \begin{subfigure}[b]{0.3\textwidth}
+        \centering
+        \includegraphics[width=\textwidth]{images/color_map_1_imag.pdf}
+        \caption{Im$(u_{\ivec, \jvec})$ at time $t=0.001$.}
+        \label{fig:colormap_1_imag}
+    \end{subfigure}
+    \hfill
+    \begin{subfigure}[b]{0.3\textwidth}
+        \centering
+        \includegraphics[width=\textwidth]{images/color_map_2_imag.pdf}
+        \caption{Im$(u_{\ivec, \jvec})$ at time $t=0.002$.}
+        \label{fig:colormap_2_imag}
+    \end{subfigure}
+    \caption{The time evolution of the probability function $p_{\ivec, \jvec}^{n}$.}
+    \label{fig:colormap}
+\end{figure*}
 \end{document}

latex/sections/conclusion.tex

@@ -2,24 +2,36 @@
 
 \begin{document}
 \section{Conclusion}\label{sec:conclusion}
+% We have simulated the two-dimensional time-dependent Schrödinger equation, to study 
+% variations of the double-slit experiment. To derive a discretized equation 
+% we applied the Crank-Nicolson scheme in 2+1 dimensions. In addition, we have used 
+% Dirichlet boundary conditions to express the equation in matrix form, and solve 
+% it using the sparse matrix solver \verb|superlu|. Our implementation, and choice 
+% of solver method, resulted in a deviation from conserved total probability on the 
+% scale $10^{-14}$, for both the single and double slit setup. To illustrate the time 
+evolution of the probability function, we created colormap plots for time steps 
+$t = [0, 0.001, 0.002]$. We also included separate plots for each time step of 
+Re$(u_{\ivec, \jvec})$ and Im$(u_{\ivec, \jvec})$. In addition, we determined the 
+normalized particle detection probability $p(y \ | \ x=0.8, t=0.002)$, for single-, 
+double- and triple-slit.
 % Rewrite this section to differ from the abstract
-We have simulated the two-dimensional time-dependent Schrödinger equation, to study 
+We simulated the two-dimensional time-dependent Schrödinger equation, and studied 
 variations of the double-slit experiment. To solve the partial differential equations 
-we have applied the Crank-Nicolson scheme in 2+1 dimensions, and derived a discretized 
+we applied the Crank-Nicolson scheme in 2+1 dimensions, and derived a discretized 
-equation. In addition, we have used Dirichlet boundary conditions to express the 
+equation. We used Dirichlet boundary conditions to simplify the equation, 
-equation in matrix form and solve it using the sparse matrix solver \verb|superlu|. 
+expressed the equation in matrix form and solved it using the sparse matrix solver 
-Our implementation, and choice of solver method, resulted in conserved total 
+\verb|superlu|. The total probability $\sum_{\ivec, \jvec} p_{\ivec, \jvec}^{n}$ 
-probability $\sum_{\ivec, \jvec} p_{\ivec, \jvec}^{n}=1$ for both the single and 
+deviated from $1.0$ by a factor of $10^{-14}$ for both the single and double slit 
-double slit setup. 
+setup. % Add something about computational accuracy?
 
-To illustrate the time evolution of the probability function $p_{\ivec, \jvec}^{n} = u_{\ivec, \jvec}^{n*} u_{\ivec, \jvec}^{n}$, 
+We illustrated the time evolution of the probability function $p_{\ivec, \jvec}^{n} = u_{\ivec, \jvec}^{n*} u_{\ivec, \jvec}^{n}$, 
-we created colormap plots at time steps $t = [0, 0.001, 0.002]$. Since we are working 
+using colormap plots for time steps $t = [0, 0.001, 0.002]$. In addition, we included 
-with complex numbers, we included separate plots for each time step of Re$(u_{\ivec, \jvec})$ 
+separate plots for each time step of Re$(u_{\ivec, \jvec})$ and Im$(u_{\ivec, \jvec})$, 
-and Im$(u_{\ivec, \jvec})$.
+to show the components of the complex values. This resulted in a visible diffraction 
-% We observed something...
+patterns for the double-slit experiment.
-
-In addition, we determined the normalized particle detection probability $p(y \ | \ x=0.8, t=0.002)$, 
-for single-, double- and triple-slit.
-% We observed something here as well...
 
+In addition, we studied the normalized particle detection probability $p(y \ | \ x=0.8, t=0.002)$, 
+for single-, double- and triple-slit. We found that increasing the number of slits 
+in the barrier, increased the number of areas of both high and low probability of 
+particle detection. 
 \end{document}

latex/sections/introduction.tex

@@ -2,6 +2,10 @@
 
 \begin{document}
 \section{Introduction}\label{sec:introduction}
+% Light: wave particle
+% Wave equation
+% Hyugens theory
+% Thomas Young
 The nature of light has long been a subject of interest and discussion. %
 % Important part of human behavior is observing and understanding our surroundings. 
 % Many big discoveries have been made through observations, verified by mathematical 

latex/sections/methods.tex

@@ -11,62 +11,64 @@ The Schrödinger equation has a general form
     i \hbar \frac{\partial}{\partial t} | \Psi \rangle &= \hat{H} | \Psi \rangle \ ,
     \label{eq:schrodinger_general}
 \end{align}
-where $i$ is the imaginary number, and $\hbar$ is Plancks constant. $\hat{H}$ is 
+where $i$ is the imaginary unit, and $\hbar$ is Plancks constant. $\hat{H}$ is 
-a Hamiltonian operator, which represent the energy for the system, and $| \Psi \rangle$ 
+a Hamiltonian operator, which represents the energy for the system, and $| \Psi \rangle$ 
 is the quantum state. In two-dimensional position space, the quantum state can 
 be expressed using the time-dependent complex-valued wave function $\Psi (x, y, t)$. 
 Using Born rule, the square modulus of the wave function is proportional to the 
-probability density of finding a particle at position $(x, y)$ at time t. The relation 
+probability density of detecting a particle at position $(x, y)$ at time $t$. The 
-is given by 
+relation is given by 
 \begin{align}
     p(x, y \ | \ t) &= |\Psi(x, y, t)|^{2} = \Psi^{*}(x, y, t) \Psi(x, y, t) \ ,
     \label{eq:born_rule}
 \end{align}
 where $\Psi^{*}$ denotes the complex conjugated wave function. 
 % Add something about kinetic and potential energy, to introduce the potential V
-When the potential is time-independent, the Schrödinger equation can be expressed as 
+When the potential is time-independent, and the particle is non-relativistic, 
+the Schrödinger equation can be expressed as 
 \begin{align*}
     i \hbar \frac{\partial}{\partial t} \Psi (x, y, t) &= - \frac{\hbar^{2}}{2m} \bigg( \frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}} \bigg) \Psi (x, y, t) \\
     & \quad + V(x, y, t) \Psi (x, y, t) \numberthis \ .
     \label{eq:schrodinger_special}
 \end{align*}
-The partial derivatives (...) gives the kinetic energy, and the potental $V$ is 
+The partial derivatives are expressions of the kinetic energy, and the potental $V$ 
-the external environment. In this experiment we will only consider the case where 
+encodes the external environment. In this experiment we will only consider the case where 
 the potential is time-independent, resulting in $V = V(x, y)$
+
 When we scale Schrödinger equation by the dimensionful variables, we are left with 
-a wave function $u$, potential $v$ and the dimensionless equation
+the wave function $u$ and the potential $v$. The dimensionless equation is given by 
 \begin{align}
     i \frac{\partial u}{\partial t} &= - \frac{\partial^{2} u}{\partial x^{2}} - \frac{\partial^{2} u}{\partial y^{2}} + v(x, y) u \ .
     \label{eq:schrodinger_dimensionless}
 \end{align} % 
-This gives us the Born rule 
+As a result of working in position space, the Born rule is given by 
 \begin{align}
-    p(x, y \ | \ t) &= |u(x, y, t)|^{2} = u^{*}(x, y, t) u(x, y, t) \ .
+    p(x, y \ | \ t) &= |u(x, y, t)|^{2} = u^{*}(x, y, t) u(x, y, t) \ ,
     \label{eq:born_rule_scaled}
 \end{align}
-
+where we assume a normalized wave function $u(x, y, t)$.
 
 
 \subsection{The Crank-Nicolson scheme}\label{ssec:crank_nicolson} %
-When we evaluate a particle in position space, we have to consider partial differential 
+When we evaluate a particles position, we have to consider partial differential 
 equations (PDE). To solve these numerically, we have to discretize Equation \eqref{eq:schrodinger_dimensionless}. 
 We use the $\theta$-rule \footnote{Using the $\theta$-rule, we can derive Forward Euler using $\theta = 1$, and Backward Euler using $\theta = 0$}, 
-to combine the forward (explicit) and backward (implicit) finite difference method. 
+to combine the forward (explicit) and backward (implicit) finite difference methods. 
 The result is a linear combination of the explicit and implicit scheme, given by 
 \begin{align}
     \frac{u_{\ivec, \jvec}^{n+1} - u_{\ivec, \jvec}^{n}}{\Delta t} &= \theta F_{\ivec, \jvec}^{n+1} + (1 - \theta) F_{\ivec, \jvec}^{n} \ ,
     \label{eq:theta_rule}
 \end{align} %
 where $\theta \in [0, 1]$.
-To simplify notation and avoid confusion of indices with the imaginary number $i$, 
+To simplify the notation, and avoid any confusion of the indices with the imaginary unit $i$, 
 we have used the notation $\ivec, \jvec$ in subscript to indicate the commonly named indices $i, j$ 
 in x- and y-direction. In addition, the superscript $n, n+1$ indicate position in time. 
-We get the Crank-Nicolson method (CN) when $\theta = 1/2$ gives Crank-Nicolson
+We derive the Crank-Nicolson scheme (CN) by using $\theta = 1/2$, given by
 \begin{align}
     \frac{u_{\ivec, \jvec}^{n+1} - u_{\ivec, \jvec}^{n}}{\Delta t} &= \frac{1}{2} \bigg[ F_{\ivec, \jvec}^{n+1} + F_{\ivec, \jvec}^{n} \bigg] \ .
-    \label{eq:crank_nicolson_method}
+    \label{eq:crank_nicolson_scheme}
 \end{align} %
-Using CN, we derive the discretized Schrödinger equation given by 
+Using CN, we derive the discretized Schrödinger equation, given by 
 \begin{align*}
     & u_{\ivec, \jvec}^{n+1} - r \big[ u_{\ivec+1, \jvec}^{n+1} - 2u_{\ivec, \jvec}^{n+1} + u_{\ivec-1, \jvec}^{n+1} \big] \\
     & - r \big[ u_{\ivec, \jvec+1}^{n+1} - 2u_{\ivec, \jvec}^{n+1} + u_{\ivec, \jvec-1}^{n+1} \big] + \frac{i \Delta t}{2} v_{\ivec, \jvec} u_{\ivec, \jvec}^{n+1} \\
@@ -76,34 +78,34 @@ Using CN, we derive the discretized Schrödinger equation given by
 \end{align*} %
 where $r$ is defined as
 \begin{align*}
-    r \equiv \frac{i \Delta t}{2 \Delta h^{2}}
+    r \equiv \frac{i \Delta t}{2 \Delta h^{2}} \ .
 \end{align*} %
-The full derivation of both Equation \eqref{eq:crank_nicolson_method} and Equation \eqref{eq:schrodinger_discretized} 
+The full derivation of both Equation \eqref{eq:crank_nicolson_scheme} and Equation \eqref{eq:schrodinger_discretized} 
 can be found in Appendix \ref{ap:crank_nicolson}.
 
 
 \subsection{The double-slit experiment}\label{ssec:double_slit} % 
 Thomas Young first performed the double-slit experiment in 1801 to demonstrate the  
 principle of interference of light \cite{britannica:2023:young}, while postulating 
-light as waves rather than particles. The double-slit experiment result in a diffraction 
+light as waves rather than particles. The double-slit experiment results in a diffraction 
 pattern on a detector screen, where constructive interference of light result in 
-bright spots, and destructive interference result in dark spots as showed in Figure 
+bright spots, and destructive interference result in dark spots. An illustration 
-\ref{fig:youngs_double_slit}.
+of Thomas Young's setup can be found in Figure \ref{fig:youngs_double_slit}.
 \begin{figure}
     \centering 
     \includegraphics[width=\linewidth]{images/youngs_double_slit.pdf}
     \caption{The setup of Thomas Young's double slit experiment, where $S_{0}$ denotes 
-    the light source, $S_{1}$ and $S_{2}$ denotes the slits in the wall.}
+    the light source, $S_{1}$ and $S_{2}$ denotes the slits in the barrier \cite[p. 4]{mit:2004:physics}.}
     \label{fig:youngs_double_slit}
 \end{figure}
 
-After the wave passes through the two slits, the pattern observed is determined by 
+After the wave passes through the barrier, the pattern observed is determined by 
-the path difference determined by 
+the path difference given by 
 \begin{align*}
     \delta = d \sin (\theta) = m \lambda \ ,
 \end{align*}
 where $\lambda$ is the wavelength and $m$ is called the order number. $d$ is the 
-distance between the center of the two slits, assuming that the distance between 
+distance between the center of the two slits, while assuming that the distance between 
 the wall and the detector screen $L >> \delta$ \cite[p. 6]{mit:2004:physics}. In 
 this case, we observe constructive interference when 
 \begin{align*}
@@ -111,24 +113,32 @@ this case, we observe constructive interference when
 \end{align*} 
 and destructive interference when
 \begin{align*}
-    \delta = (m + \frac{1}{2}) \lambda && m = 0, \pm 1, \pm 2 \dots \ ,
+    \delta = (m + \frac{1}{2}) \lambda && m = 0, \pm 1, \pm 2 \dots \ .
 \end{align*} 
-
 % Something about Heisenberg uncertainty principle 
 
+
 \subsection{Implementation}\label{ssec:implementation} % 
-% Add tables of parameters used, initial conditions, notation etc.
+In this experiment, we set up the grid with an equal step size in x- and y-direction $h$, 
-A, B are sparse csc matrix
+and step size in t-direction $\Delta t$, such that
-- theory of csc matrix?
+\begin{align*}
-% \begin{equation*}
+    x \in [0, 1] && x \rightarrow x_{\ivec} = \ivec h && \ivec = 0, 1, \dots, M-1 \\
-% \begin{pNiceArray}{ccc|ccc}
+    y \in [0, 1] && y \rightarrow y_{\jvec} = \jvec h && \jvec = 0, 1, \dots, M-1 \\
-%     \bullet & \bullet &  & \bullet &  &  &  &  &  \\
+    t \in [0, T] && t \rightarrow t_{n} = n \Delta t && n = 0, 1, \dots, N_{t}-1 \ .
-%     \bullet & \bullet & \bullet &  & \bullet &  &  &  & \\
+\end{align*}
-%       & \bullet & \bullet &  &  & \bullet &  &  & 
+In addition, we simplified the indices such that 
-% \end{pNiceArray}
+\begin{align*}
-% \end{equation*}
+    u(x, y, t) \rightarrow u(\ivec h, \jvec h, n \Delta t) \equiv u_{\ivec, \jvec}^{n} \\
-We use Dirichlet boundary conditions, as given in Table \ref{tab:boundary_conditions}, 
+    v(x, y) \rightarrow u(\ivec h, \jvec h) \equiv v_{\ivec, \jvec} \ ,
-which allows us to express Equation \eqref{eq:schrodinger_discretized} as a matrix 
+\end{align*}
+which results in a matrix $U^{n}$ that contains elements $u_{\ivec, \jvec}^{n}$, and 
+a matrix $V$ that contains elements $v_{\ivec, \jvec}$. We used Dirichlet boundary 
+conditions, given by 
+\begin{align*}
+    u(x=0, y, t) &= 0  & u(x=1, y, t) &= 0 \\
+    u(x, y=0, t) &= 0 & u(x, y=1, t) &= 0 \ ,
+\end{align*}
+which allowed us to express Equation \eqref{eq:schrodinger_discretized} as a matrix 
 equation 
 \begin{align}
     A u^{n+1} = B u^{n} \ .
@@ -137,8 +147,8 @@ Here, both $u^{n+1}$ and $u^{n}$ are column vectors containing the internal poin
 of the $xy$ grid at time step $n+1$ and $n$, respectively. Since we have $M$ points 
 in $x$- and $y$-direction, we have $M-2$ internal points. Both $u$ vectors have 
 length $(M-2)^{2}$, and the matrices $A$ and $B$ have size $(M-2)^{2} \times (M-2)^{2}$. 
-The matrices can be decomposed as submatrices of size $(M-2) \times (M-2)$, with 
+The matrices are sparse and can be decomposed as submatrices of size $(M-2) \times (M-2)$, 
-the following pattern 
+with the following pattern 
 \begin{align*}
     A, B = 
     \begin{bmatrix}
@@ -195,7 +205,7 @@ the following pattern
            \bullet & \bullet & \bullet \\
            \phantom{\bullet} & \bullet & \bullet 
        \end{matrix}
-    \end{bmatrix}
+    \end{bmatrix} \ .
 \end{align*}
 To fill the matrices $A$ and $B$, we used 
 \begin{align*}
@@ -204,37 +214,16 @@ To fill the matrices $A$ and $B$, we used
 \end{align*} 
 An example of filled matrices can be found in Appendix \ref{ap:matrix_structure}.
 
-Notations:
+For the general setup of the barrier, we used the values in Table \ref{tab:barrier_setup}, 
-In addition, we use an equal step size in x- and y-direction, $h$ such that 
+and for the simulations, we used the parameter settings in Table \ref{tab:sim_settings}. 
+In addition, we initialized the wave function using a Gaussian wavepacket, given by 
 \begin{align*}
-    x \in [0, 1] && x \rightarrow x_{\ivec} = \ivec h && \ivec = 0, 1, \dots, M-1 \\
+    u(x, y, t=0) &= e^{- \frac{(x-x_{c})^{2}}{2 \sigma_{x}^{2}} - \frac{(y-y_{c})^{2}}{2 \sigma_{y}^{2}} + ip_{x}x + ip_{y}y} \ .
-    y \in [0, 1] && y \rightarrow y_{\jvec} = \jvec h && \jvec = 0, 1, \dots, M-1 \\
-    t \in [0, T] && t \rightarrow t_{n} = n \Delta t && n = 0, 1, \dots, N_{t}-1 
 \end{align*}
-And simplify indices such that 
+$x_{c}$ and $y_{c}$ are the coordinates of the center of the wavepacket, $\sigma_{x}$ 
-\begin{align*}
+and $\sigma_{y}$ are the width of the wavepacket, when it is initialized. The wave 
-    u(x, y, t) \rightarrow u(\ivec h, \jvec h, n \Delta t) \equiv u_{\ivec, \jvec}^{n} \\
+packet momenta are given by $p_{x}$ and $p_{y}$. 
-    v(x, y) \rightarrow u(\ivec h, \jvec h) \equiv v_{\ivec, \jvec}
+% Insert Heisenberg uncertainty here? Or refer to it?
-\end{align*}
-which gives a matrix $U^{n}$ that contains elements $u_{\ivec, \jvec}^{n}$, and 
-a matrix $V$ that contains elements $v_{\ivec, \jvec}$.
-\begin{table}[H]
-    \centering
-    \begin{tabular}{l r} % @{\extracolsep{\fill}}
-        \hline
-        Position & Value  \\
-        \hline 
-        $u(x=0, y, t)$ & $0$  \\
-        $u(x=1, y, t)$ & $0$ \\
-        $u(x, y=0, t)$ & $0$  \\
-        $u(x, y=1, t)$ & $0$  \\
-        \hline
-    \end{tabular}
-    \caption{Boundary conditions in the xy-plane, also known as Dirichlet boundary conditions.}
-    \label{tab:boundary_conditions}
-\end{table}
-
-For the general setup of the wall, we used 
 \begin{table}[H]
     \centering
     \begin{tabular}{l r} % @{\extracolsep{\fill}}
@@ -247,15 +236,14 @@ For the general setup of the wall, we used
         Slit aperture & $0.05$  \\
         \hline
     \end{tabular}
-    \caption{Wall setup.}
+    \caption{Barrier parameters and values.}
-    \label{tab:wall_setup}
+    \label{tab:barrier_setup}
 \end{table}
-
 \begin{table}[H]
     \centering
     \begin{tabular}{l r r} % @{\extracolsep{\fill}}
         \hline
-        Simulation & $1$ & $2$ \\
+        Parameter & Setting 1 & Setting 2 \\
         \hline 
         $h$ & $0.005$ & $0.005$ \\
         $\Delta t$ & $2.5 \times 10^{-5}$ & $2.5 \times 10^{-5}$ \\
@@ -269,12 +257,39 @@ For the general setup of the wall, we used
         $v_{0}$ & $0$ & $1 \times10^{10}$ \\
         \hline
     \end{tabular}
-    \caption{Wall setup.}
+    \caption{Simulation settings used in the double slit experiment. Setting 1 is 
-    \label{tab:sim_setup}
+    used when the barrier is switched off and setting 2 is used when the barrier 
+    switched on.}
+    \label{tab:sim_settings}
 \end{table}
 
+To check if the total probability is conserved over time, and that the implementation 
+was correct, we computed the deviation from $1.0$ given by
+\begin{align*}
+    s^{n} &= |\sum_{\ivec , \jvec} p_{\ivec , \jvec}^{n} - 1| \\
+    &= |\sum_{\ivec , \jvec} u_{\ivec , \jvec}^{n*} u_{\ivec , \jvec}^{n} - 1| \ .
+\end{align*}
+
 \subsection{Tools}\label{ssec:tools} %
 The double-slit experiment is implemented in C++. We use the Python library 
 \verb|matplotlib| \cite{hunter:2007:matplotlib} to produce all the plots, and 
 \verb|seaborn| \cite{waskom:2021:seaborn} to set the theme in the figures. 
 \end{document}
+
+
+
+% \begin{table}[H]
+%     \centering
+%     \begin{tabular}{l r} % @{\extracolsep{\fill}}
+%         \hline
+%         Position & Value  \\
+%         \hline 
+%         $u(x=0, y, t)$ & $0$  \\
+%         $u(x=1, y, t)$ & $0$ \\
+%         $u(x, y=0, t)$ & $0$  \\
+%         $u(x, y=1, t)$ & $0$  \\
+%         \hline
+%     \end{tabular}
+%     \caption{Boundary conditions in the xy-plane, also known as Dirichlet boundary conditions.}
+%     \label{tab:boundary_conditions}
+% \end{table}

latex/sections/results.tex

@@ -4,32 +4,89 @@
 \section{Results}\label{sec:results}
 \subsection{Deviation}\label{ssec:deviation}
 % Problem 3: Discuss approaches to solve Au^{n+1} = b, dealing with sparse matrix...
-We used the superlu solver, which is a dedicated solver for sparse matrices. It is 
+We used the \verb|superlu| solver, which is a solver for sparse matrices. It is 
-generally used to solve nonsymmetric, sparse matrices. However, as the lapack solver 
+generally used to solve nonsymmetric, sparse matrices. However, as the alternative 
-converts the sparse matrix to a dense matrix, it will increase memory usage compared 
+solver \verb|lapack| converts a sparse matrix to a dense matrix, it will increase 
-to superlu.
+memory usage compared to \verb|superlu|.
 % Problem 7: Consequenses of solver choice, in regards to accuracy of probability conserved
 %            Add plot of deviation for both single- and double-slit
-Since we use a solver for sparse matrices, we decrease number of computations performed 
+Since we used a solver for sparse matrices, we decrease the number of computations performed
-compared to solver using dense matrix. We check if the total probability is conserved 
+compared to number of computations using a solver for dense matrices. 
-over time, by plotting the deviation $s$ as
+We checked if the total probability was conserved over time, by plotting the deviation 
-\begin{align*}
+from $1.0$.
-    s^{n} = 1 - \sum_{\ivec , \jvec} p_{\ivec , \jvec}^{n} = 1 - \sum_{\ivec , \jvec} u_{\ivec , \jvec}^{n*} u_{\ivec , \jvec}^{n} \ .
-\end{align*}
-The deviation as a function of time is plotted in Figure \ref{fig:deviation}.
 \begin{figure}
     \centering
     \includegraphics[width=\linewidth]{images/probability_deviation.pdf}
-    \caption{Deviation for $t \in [0, T]$ where $T=0.008$.}
+    \caption{Deviation of total probability, for time $t \in [0, T]$ where $T=0.008$.}
     \label{fig:deviation}
 \end{figure}
+In Figure \ref{fig:deviation}, we observe a larger deviation of total probability 
+for a barrier with double slits compared to no barrier. Interaction with the 
+barrier result in a change ... of values, as the light passes through, 
+The result is more prone to computational errors. Using no barrier, the light is not 
+affected by obstacles resulting in a more stable deviation from the total probability. 
+In addition, we have to consider the limitation of a computer, some computational 
+error is to be expected.
 
 \subsection{Time evolution}\label{ssec:time_evolution}
 % Problem 8: Colormap, include plot of both Re and Im for different time steps
 %            Account for color scale
+We ran the simulation using the values in column 2 of Table \ref{tab:sim_setup}, found 
+in Section \ref{ssec:implementation}. To study the time evolution of the probability 
+function, we created colormap plots for different time steps. Figure \ref{fig:colormap_0_prob}, 
+Figure \ref{fig:colormap_1_prob}, and Figure \ref{fig:colormap_2_prob} show the 
+result for time steps $t=[0, 0.001, 0.002]$, respectively. In addition, we created 
+separate plots for the real and imaginary part of $u_{\ivec, \jvec}$, for the same 
+time steps. The result can be found in Appendix \ref{ap:figures}, in Figure \ref{fig:colormap}.
+\begin{figure}
+    \centering
+    \includegraphics[width=\linewidth]{images/color_map_0_prob.pdf}
+    \caption{The probability function $p_{\ivec, \jvec}^{n}$, at time $t=0$.}
+    \label{fig:colormap_0_prob}
+\end{figure}
+\begin{figure}
+    \centering
+    \includegraphics[width=\linewidth]{images/color_map_1_prob.pdf}
+    \caption{The probability function $p_{\ivec, \jvec}^{n}$, at time $t=0.001$.}
+    \label{fig:colormap_1_prob}
+\end{figure}
+\begin{figure}
+    \centering
+    \includegraphics[width=\linewidth]{images/color_map_2_prob.pdf}
+    \caption{The probability function $p_{\ivec, \jvec}^{n}$, at time $t=0.002$.}
+    \label{fig:colormap_2_prob}
+\end{figure}
+In Figure \ref{fig:colormap_1_prob}, the ... interacts with the double 
+slitted barrier, which result in a wavelike pattern. However, the waves are more 
+visible when we observe the real and imaginary part separately in Figure \ref{fig:colormap}.
+
 
 \subsection{Particle detection}\label{ssec:particle_detection}
 % Problem 9: Plot detection probability for single-, double- and triple-slit
-
+We used the simulation from the previous section, assumed a detector screen 
+located at $x=0.8$, and plotted the detection probability along the screen at time 
+$t=0.002$. We adjusted the parameters to include single-, double-, and triple-slitted 
+barrier. The results is found in Figure \ref{}, Figure \ref{}, and Figure \ref{fig:}.
+\begin{figure}
+    \centering
+    \includegraphics[width=\linewidth]{images/single_slit_detector.pdf}
+    \caption{Probability of particle detection along a detector screen at time $t=0.002$, 
+    when using a single-slit barrier.}
+    \label{fig:particle_detection_single}
+\end{figure}
+\begin{figure}
+    \centering
+    \includegraphics[width=\linewidth]{images/double_slit_detector.pdf}
+    \caption{Probability of particle detection along a detector screen at time $t=0.002$, 
+    when using a double-slit barrier.}
+    \label{fig:particle_detection_double}
+\end{figure}
+\begin{figure}
+    \centering
+    \includegraphics[width=\linewidth]{images/triple_slit_detector.pdf}
+    \caption{Probability of particle detection along a detector screen at time $t=0.002$, 
+    when using a triple-slit barrier.}
+    \label{fig:particle_detection_triple}
+\end{figure}
 
 \end{document}

python_scripts/normalization.py

@@ -21,9 +21,17 @@ def plot():
         "data/probability_deviation_no_slits.txt",
         "data/probability_deviation_slits.txt",
     ]
-    for file in files:
+    labels = [
+        "none",
+        "double-slit"
+    ]
+    colors = sns.color_palette("mako", n_colors=2)
+
+    fig, ax = plt.subplots() 
+
+    for file, label, color in zip(files, labels, colors):
         with open(file) as f:
-            lines = f.readlines();
+            lines = f.readlines()
             x = []
             arr = []
             for line in lines:
@@ -31,10 +39,13 @@ def plot():
                 x.append(float(tmp[0]))
                 arr.append(float(tmp[1]))
 
-            plt.plot(x,arr)
+            ax.plot(x, arr, label=label, color=color)
+        
+        ax.legend(title="Barrier")
+        ax.set_xlabel("Time (t)")
+        ax.set_ylabel("Deviation")
+    fig.savefig("latex/images/probability_deviation.pdf", bbox_inches="tight")
 
-        plt.savefig("latex/images/probability_deviation.pdf")
 
 if __name__ == "__main__":
-    plot()
+    plot()
-