Merge pull request #16 from FYS3150-G2-2023/develop

Develop

.gitignore

@@ -37,6 +37,7 @@
 *.out
 *.synctex.gz
 *.bbl 
+*.blg
 latex/mainNotes.bib
 
 

latex/appendix.tex

@@ -0,0 +1,266 @@
+\documentclass[../main.tex]{subfiles}
+\graphicspath{{\subfix{../images/}}}
+
+\begin{document}
+\section{Derivation of equations}\label{sec:derivations}
+% Problem 1
+\subsection{Equations of motion}\label{sec:eq_motion} %
+First, we need to define the velocity of the particle
+\begin{align*}
+    \mathbf{v} \equiv \frac{d \mathbf{r}}{dt} &= \bigg( \frac{dx}{dt}, \frac{dy}{dt}, \frac{dz}{dt} \bigg).
+\end{align*} %
+%
+We can rewrite the velocity as $\dot{r} = (\dot{x}, \dot{y}, \dot{z})$, and find the cross product
+\begin{align*}
+    q \mathbf{v} \cross \mathbf{B} &= 
+    q \begin{vmatrix}
+        \hat{e}_{x} & \hat{e}_{y} & \hat{e}_{z} \\
+        \dot{x} & \dot{y} & \dot{z} \\
+        0 & 0 & B_{0}
+    \end{vmatrix}
+    = q \big( B_{0} \dot{y}, -B_{0} \dot{x}, 0 \big).
+\end{align*} %
+%
+We are considering an ideal Penning traps, where we define the electric potential as
+\begin{align*}
+    V(x, y, z) &= \frac{V_{0}}{2 d^{2}}(2z^{2} - x^{2} - y^{2}).
+\end{align*} %
+%
+The relationship between the electric field $\mathbf{E}$ and the electric potential of the field is given by
+\begin{align*}
+    \mathbf{E} &= - \nabla V \\
+    &= - \bigg( \frac{dV}{dx}, \frac{dV}{dy} \frac{dV}{dz} \bigg) \\
+    &= \frac{V_{0}}{d^{2}} \big( x, y, -2z \big).
+\end{align*} %
+%
+We can now express the Lorentz force as
+\begin{align*}
+    \mathbf{F} &= q \mathbf{E} + q \mathbf{v} \cross \mathbf{B} \\
+    &= \frac{q V_{0}}{d^{2}} \big( x, y, -2z \big) + \big(q B_{0} \dot{y}, -q B_{0} \dot{x}, 0 \big),
+\end{align*} %
+%
+and insert it into Newtons equation \eqref{eq:newton_second}. We get
+\begin{align*}
+    \ddot{\mathbf{r}} &= \bigg( \frac{q V_{0}}{m d^{2}} x, \frac{q V_{0}}{m d^{2}} y, -\frac{2 q V_{0}}{m d^{2}} z \bigg) + \bigg(\frac{q B_{0}}{m} \dot{y}, -\frac{q B_{0}}{m} \dot{x}, 0 \bigg),
+\end{align*} %
+%
+which can be written as
+\begin{align*}
+    \ddot{x} &= \frac{q V_{0}}{m d^{2}} x + \frac{q B_{0}}{m} \dot{y}, \\
+    \ddot{y} &= \frac{q V_{0}}{m d^{2}} y - \frac{q B_{0}}{m} \dot{x}, \\
+    \ddot{z} &= -\frac{2 q V_{0}}{m d^{2}} z.
+\end{align*} %
+%
+If we define
+\begin{equation*}
+    \omega_{0} \equiv \frac{q B_{0}}{m}, \quad \omega_{z}^{2} \equiv \frac{2 q V_{0}}{m d^{2}},
+\end{equation*} %
+%
+the equations of motion can be written as
+\begin{align*}
+    \ddot{x} &= \frac{1}{2} \omega_{z}^{2} x + \omega_{0} \dot{y}, \\
+    \ddot{y} &= \frac{1}{2} \omega_{z}^{2} y - \omega_{0} \dot{x}, \\
+    \ddot{z} &= -\omega_{z}^{2} z. \\
+\end{align*} %
+%
+\subsection{General solution}\label{sec:eq_general}
+We consider the characteristic equation of a second order differential equation \cite{lindstrom:2016:ch10:5},  
+\begin{align*}
+    r^{2} + \omega_{z}^{2} &= 0 \\
+    r &= \pm \sqrt{- \omega_{z}^{2}}.
+\end{align*} %
+The characteristic equation has two complex roots
+\begin{equation*}
+     r_{1} = - i \omega_{z}, \quad r_{2} = i \omega_{z},
+\end{equation*}
+which give us solutions in the general form
+\begin{equation*}
+    z = c_{1} e^{i \omega_{z} t} + c_{2} e^{-i \omega_{z} t}.
+\end{equation*}
+In addition, for a complex number $z = a + ib$, we can define $e^{z} \equiv e^{a} (\cos{b} + i \sin{b})$ \cite{lindstrom:2016:ch3}. We can rewrite the general solution as
+\begin{align*}
+    c_{1} e^{i \omega_{z} t} + c_{2} e^{-i \omega_{z} t} &= c_{1} (\cos{\omega_{z} t} + i \sin{\omega_{z} t}) \\
+    & \quad + c_{2} (\cos{\omega_{z} t} - i \sin{\omega_{z} t} \\
+    &= E \cos{\omega_{z} t} + i F \sin{\omega_{z} t}.
+\end{align*} %
+%
+\subsection{Complex function}\label{sec:eq_complex} %
+In sec. \ref{sec:eq_motion} we found the differential equations for $\ddot{x}$ and $\ddot{y}$. To derive a single differential equation, we introduce the complex function $f(t) = x(t) + iy(t)$, which gives us
+\begin{align*}
+    0 &= \Big(\ddot{x} - \omega_{0} \dot{y} - \frac{1}{2} \omega_{z}^{2} x \Big) + i \Big(\ddot{y} + \omega_{0} \dot{x} - \frac{1}{2} \omega_{z}^{2} y \Big) \\
+    &= \ddot{x} - \omega_{0} \dot{y} - \frac{1}{2} \omega_{z}^{2} x + i\ddot{y} + i\omega_{0} \dot{x} - i \frac{1}{2} \omega_{z}^{2} y \\
+    &= \ddot{x} + i\ddot{y} + i\omega_{0} \dot{x} - \omega_{0} \dot{y} - \frac{1}{2} \omega_{z}^{2} x - i \frac{1}{2} \omega_{z}^{2} y. \\ 
+\end{align*}
+Using the definition $i = \sqrt{-1}$, we can rewrite 
+\begin{equation*}
+    i \omega_{0} \dot{x} + (-1) \omega_{0} \dot{y} = i \omega_{0} \dot{x} + i^{2} \omega_{0} \dot{y}.
+\end{equation*} 
+This gives us a single differential equation
+\begin{align*}
+    0 &= \ddot{x} + i\ddot{y} + i\omega_{0} (\dot{x} + i \dot{y}) - \frac{1}{2} \omega_{z}^{2} x - i \frac{1}{2} \omega_{z}^{2} y \\
+    &= \ddot{f} + i \omega_{0} \dot{f} - \frac{1}{2} \omega_{z}^{2} f.
+\end{align*}
+%
+\subsection{Physical coordinates}\label{sec:eq_coord}
+We can rewrite eq. \eqref{eq:general_solution}, as
+\begin{align*}
+    f(t) &= A_{+}e^{-i(\omega_{+} t + \phi_{+})} + A_{-}e^{-i(\omega_{-} t + \phi_{-})} \\
+    &= A_{+}(\cos{(\omega_{+} t + \phi_{+})} - i \sin{(\omega_{+} t + \phi_{+})}) \\
+    % \numberthis \label{eq:general_solution_trig}
+    & \quad + A_{-}(\cos{(\omega_{-} t + \phi_{-})} - i \sin{(\omega_{-} t + \phi_{-})}).
+\end{align*} %
+%
+\subsection{Upper and lower bounds}\label{sec:upper_lower_bound}
+To obtain the upper and lower bounds of the particle's distance from the origin, we first find an expression for the second norm defined as $|f(t)| = \sqrt{(x(t))^{2} + (y(t))^{2}}$. 
+\begin{align*}
+    (x(t))^{2} &= \big( A_{+}\cos(\omega_{+} t + \phi_{+}) + A_{-}\cos(\omega_{-} t + \phi_{-}) \big)^{2} \\
+    &= A_{+}^{2} \cos^{2}(\omega_{+} t + \phi_{+}) \\
+    & \quad + 2 A_{+}A_{-} \cos(\omega_{+} t + \phi_{+})\cos(\omega_{-} t + \phi_{-}) \\
+    & \quad + A_{-}^{2}\cos^{2}(\omega_{-} t + \phi_{-}), \\
+\end{align*} %
+\begin{align*}
+    (y(t))^{2} &= \big( - A_{+} \sin(\omega_{+} t + \phi_{+}) - A_{-} \sin(\omega_{-} t + \phi_{-}) \big)^{2} \\
+    &= A_{+}^{2} \sin^{2}(\omega_{+} t + \phi_{+}) \\
+    & \quad + 2 A_{+}A_{-} \sin(\omega_{+} t + \phi_{+})\sin(\omega_{-} t + \phi_{-}) \\
+    & \quad + A_{-}^{2} \sin^{2}(\omega_{-} t + \phi_{-}).
+\end{align*} %
+We insert these expressions, and find 
+\begin{align*}
+    |f(t)| &= \sqrt{(x(t))^{2} + (y(t))^{2}} \\
+    &= \sqrt{A_{+}^{2} + 2 A_{+} A_{-} \cos^{2}(\alpha) + A_{-}^{2}},
+\end{align*} %
+where $\alpha = (\omega_{+} - \omega_{-}) t +( \phi_{+} -  \phi_{-})$. If we set $\alpha = 0$ we get $\cos(0) = 1$, and obtain the upper bound
+\begin{align*}
+    R_{+} &= \sqrt{A_{+}^{2} + 2 A_{+} A_{-} + A_{-}^{2}} \\
+    &= \sqrt{(A_{+} + A_{-})^{2}} \\
+    &= A_{+} + A_{-}.
+\end{align*}
+If $\alpha = \pi$ we get $\cos(\pi) = -1$, and find the lower bound
+\begin{align*}
+    R_{-} &= \sqrt{A_{+}^{2} - 2 A_{+} A_{-} + A_{-}^{2}} \\
+    &= \sqrt{(A_{+} - A_{-})^{2}} \\
+    &= |A_{+} - A_{-}|.
+\end{align*} %
+%
+\subsection{Bounded solution}\label{sec:bounded_solution}
+To find a bounded solution, we need to consider the angular rate in eq. \eqref{eq:angular_rate}. Specifically the square-root expression
+\begin{align*}
+    \sqrt{\omega_{0}^{2} - 2 \omega_{z}^{2}}.
+\end{align*}
+A bounded solution can be found when this expression is greater than zero. Meaning
+\begin{align*}
+    \omega_{0}^{2} &> 2 \omega_{z}^{2}.
+\end{align*} %
+%
+We can now use the definition of $\omega_{0}$ and $\omega_{z}^{2}$ to find an expression related to the Penning trap parameters, and the particle properties, to get
+\begin{align*}
+    \Big( \frac{q B_{0}}{m} \Big)^{2} &> 2 \frac{2 q V_{0}}{m d^{2}} \\
+    \frac{q^{2} B_{0}^{2}}{m^{2}} &> \frac{4 q V_{0}}{m d^{2}}.
+\end{align*} %
+
+\section{Values used in simulation}
+\subsection{Specific analytical solution}\label{sec:spec_analytical}
+To compare our implementation of the Penning trap and particle, we have to specify initial conditions of the system to compare with the analytical solution. The initial conditions of a particle with a single positive charge, can be found in table \ref{tab:initial_particle_cond}. 
+\begin{table}[H]
+    \centering
+    \begin{tabular}[c]{ll}
+        $\mathbf{r}(0)$ & $\dot{\mathbf{r}}(0)$ \\
+        \hline
+        $x(0)$ = $x_{0}$ & $\dot{x}(0)$ = $0$ \\
+        $y(0)$ = $0$ & $\dot{y}(0)$ = $v_{0}$ \\
+        $z(0)$ = $z_{0}$ & $\dot{z}(0)$ = $0$ \\
+        \hline
+    \end{tabular}
+    \caption{Initial values of a particle with a single charge $q$, and mass $m$, confined in a Penning trap.}
+    \label{tab:initial_particle_cond}
+\end{table}
+
+\begin{table}[H]
+    \centering
+        \begin{tabular}[c]{lll}
+            Property & Value \\
+            \hline
+            $B_{0}$ & $1.00 T$ = $9.65 \cross 10^{1} \frac{u}{(\mu s) e}$ \\
+            $V_{0}$ & $25.0 mV$ = $2.41 \cross 10^{6} \frac{u (\mu m)^{2}}{(\mu s)^{2} e}$ \\
+            $d$ & $500 \ \mu m$ = \\
+            \hline
+        \end{tabular}
+    \caption{Default configuration of the Penning trap, where the value of T and V can be found in table \ref{tab:constants}.}
+    \label{tab:penning_config}
+\end{table}
+
+\begin{table}[H]
+    \centering
+        \begin{tabular}[c]{lll}
+            Property & Value \\
+            \hline
+            $q$ & $1.00 T$ & $9.65 \cross 10^{1} \frac{u}{(\mu s) e}$ \\
+            $m$ & $25.0 mV$ & $2.41 \cross 10^{6} \frac{u (\mu m)^{2}}{(\mu s)^{2} e}$ \\
+            \hline
+        \end{tabular}
+    \caption{Default configuration of the Penning trap, where the value of T and V can be found in table \ref{tab:constants}.}
+    \label{tab:particle_config}
+\end{table}
+
+\begin{table}[H]
+    \centering
+    \begin{tabular}[c]{lll}
+        $n_{k}$ & Time steps & Step size \\
+        \hline
+        $n_{1}$ & $4000$ & $0.0125$ \\
+        $n_{2}$ & $8000$ & $0.00625$ \\
+        $n_{3}$ & $16000$ & $0.003125$ \\
+        $n_{4}$ & $32000$ & $0.0015625$ \\
+        \hline
+    \end{tabular}
+    \caption{Number of steps used in a given simulation $k$, where the step size corresponds to $h_{1} = 50 / n_{k} \mu s$.}
+    \label{tab:time_steps}
+\end{table}
+
+\subsection{Numbers and units}\label{sec:numbers_units}
+\begin{table}[H]
+    \centering
+    \begin{tabular}[c]{lll}
+        Constant & Value & Unit \\
+        \hline
+        $k_{e}$ (Coulomb) & $1.38935333 \cross 10^{5}$ & $\frac{u (\mu m)^{3}}{(\mu s)^{2} e^{2}}$ \\
+        T (Tesla) & $9.64852558 \cross 10^{1}$ & $\frac{u}{(\mu s) e}$ \\
+        V (Volt) & $9.64852558 \cross 10^{7}$ & $\frac{u (\mu m)^{2}}{(\mu s)^{2} e}$ \\
+        \hline
+    \end{tabular}
+    \caption{Value of the Coulomb constant ($k_{e}$), and the SI units for magnetic field strength ($T$) and electric potential ($V$). The base units are given by length in micrometre ($\mu m$), time in microseconds ($\mu s$), mass in ($u$), and charge in elementary charge ($e$).}
+    \label{tab:constants}
+\end{table}
+
+\section{Algorithm implementation and testing}\label{sec:algo}
+\subsection{Forward Euler}\label{sec:algo_euler}
+For a particle $i$, at time step $j$, the forward Euler method for a coupled system can be expressed as
+\begin{align*}
+    \mathbf{r}_{i,j+1} &= \mathbf{r}_{i,j} + h \frac{d \mathbf{r}_{i,j}}{dt} = \mathbf{r}_{i,j} + h \mathbf{v}_{i,j} \\
+    \mathbf{v}_{i,j+1} &= \mathbf{v}_{i,j} + h \frac{\mathbf{v}_{i,j}}{dt} = \mathbf{v}_{i,j} + h \frac{\mathbf{F}(t_{j},\mathbf{v}_{i,j}, \mathbf{r}_{i,j} )}{m_{i}},
+\end{align*}
+$m_{i}$ is the mass of the particle, and $h$ is the step length.
+
+\subsection{4th order Runge-Kutta}\label{sec:algo_rk4}
+For a particle $i$, at time step $j$, the 4th order Runge-Kutta method for a coupled system can be expressed as
+\begin{align*}
+    \mathbf{v}_{i,j+1} &= \mathbf{v}_{i,j} + \frac{h}{6} (\mathbf{k}_{\mathbf{v},1,i} + 2 \mathbf{k}_{\mathbf{v},2,i} + 2 \mathbf{k}_{\mathbf{v},3,i} + \mathbf{k}_{\mathbf{v},4,i} ), \\
+    \mathbf{r}_{i,j+1} &= \mathbf{r}_{i,j} + \frac{h}{6} (\mathbf{k}_{\mathbf{r},1,i} + 2 \mathbf{k}_{\mathbf{r},2,i} + 2\mathbf{k}_{\mathbf{r},3,i} + \mathbf{k}_{\mathbf{r},4,i}),
+\end{align*} 
+where
+\begin{align*}
+    \mathbf{k}_{\mathbf{v},1,i} &= \frac{\mathbf{F}_{i}(t_{j}, \mathbf{v}_{i,j}, \mathbf{r}_{i,j})}{m_{i}}, \\
+    \mathbf{k}_{\mathbf{r},1,i} &= \mathbf{v}_{i,j}, \\
+    \mathbf{k}_{\mathbf{v},2,i} &= \frac{\mathbf{F}_{i}(t_{j}+\frac{h}{2}, \mathbf{v}_{i,j} + h \frac{\mathbf{k}_{\mathbf{v},1,i}}{2}, \mathbf{r}_{i,j} + h \frac{\mathbf{k}_{\mathbf{r},1,i}}{2})}{m_{i}}, \\
+    \mathbf{k}_{\mathbf{r},2,i} &= \mathbf{v}_{i,j} + h \frac{\mathbf{k}_{\mathbf{v},1,i}}{2}, \\
+    \mathbf{k}_{\mathbf{v},3,i} &= \frac{\mathbf{F}_{i}(t_{j}+\frac{h}{2}, \mathbf{v}_{i,j} + h \frac{\mathbf{k}_{\mathbf{v},2,i}}{2}, \mathbf{r}_{i,j} + h \frac{\mathbf{k}_{\mathbf{r},2,i}}{2})}{m_{i}}, \\
+    \mathbf{k}_{\mathbf{r},3,i} &= \mathbf{v}_{i,j} + h \frac{\mathbf{k}_{\mathbf{v},2,i}}{2}, \\
+    \mathbf{k}_{\mathbf{v},4,i} &= \frac{\mathbf{F}_{i}(t_{j}+h, \mathbf{v}_{i,j} + h \mathbf{k}_{\mathbf{v},3,i}, \mathbf{r}_{i,j} + h \mathbf{k}_{\mathbf{r},3,i})}{m_{i}} \\
+    \mathbf{k}_{\mathbf{r},4,i} &= \mathbf{v}_{i,j} + h \frac{\mathbf{k}_{\mathbf{v},1,i}}{2}.
+\end{align*}
+In order to find each $\mathbf{k}_{\mathbf{r},i}$ and $\mathbf{k}_{\mathbf{v},i}$, we need to first compute all $\mathbf{k}_{\mathbf{r},i}$ and $\mathbf{k}_{\mathbf{v},i}$ for all particles, then update the particle values in order to compute $\mathbf{k}_{\mathbf{r},i+1}$ and $\mathbf{k}_{\mathbf{v},i+1}$. 
+
+\subsection{Test of implementation}\label{sec:testing}
+We implemented a test suite, to validate the implementation during code development. In addition, we have implemented functionality to get informative output while testing the code. Further instructions with code can be found in the github repo \cite{github:repo}.
+
+\end{document}

latex/main.blg

@@ -1,53 +0,0 @@
-This is BibTeX, Version 0.99d (TeX Live 2022/CVE-2023-32700 patched)
-Capacity: max_strings=200000, hash_size=200000, hash_prime=170003
-The top-level auxiliary file: main.aux
-The style file: unsrt.bst
-Database file #1: mainNotes.bib
-Database file #2: references/references.bib
-Warning--entry type for "britannica:2023:matter" isn't style-file defined
---line 63 of file references/references.bib
-Warning--entry type for "openmp:2018" isn't style-file defined
---line 92 of file references/references.bib
-Warning--empty chapter and pages in vogel:2018:ch1
-You've used 7 entries,
-            1791 wiz_defined-function locations,
-            485 strings with 4118 characters,
-and the built_in function-call counts, 1494 in all, are:
-= -- 143
-> -- 28
-< -- 2
-+ -- 14
-- -- 7
-* -- 72
-:= -- 212
-add.period$ -- 19
-call.type$ -- 7
-change.case$ -- 7
-chr.to.int$ -- 0
-cite$ -- 8
-duplicate$ -- 79
-empty$ -- 172
-format.name$ -- 7
-if$ -- 355
-int.to.chr$ -- 0
-int.to.str$ -- 7
-missing$ -- 8
-newline$ -- 36
-num.names$ -- 7
-pop$ -- 35
-preamble$ -- 1
-purify$ -- 0
-quote$ -- 0
-skip$ -- 46
-stack$ -- 0
-substring$ -- 95
-swap$ -- 25
-text.length$ -- 2
-text.prefix$ -- 0
-top$ -- 0
-type$ -- 0
-warning$ -- 1
-while$ -- 16
-width$ -- 8
-write$ -- 75
-(There were 3 warnings)

latex/main.tex

@@ -49,7 +49,7 @@
 
 \begin{document}
 
-\title{Project 3}  % self-explanatory
+\title{Simulating Particles in an Ideal Penning Trap}
 \author{Cory Alexander Balaton \& Janita Ovidie Sandtrøen Willumsen \\ \faGithub \, \url{https://github.uio.no/FYS3150-G2-2023/Project-3}} % self-explanatory
 \date{\today}                             % self-explanatory
 \noaffiliation                            % ignore this, but keep it.
@@ -73,11 +73,13 @@
 % Conclusion
 \subfile{sections/conclusion}
 
+\newpage
 
 \appendix
-\subfile{appendix/appendix_a}
+\subfile{appendix}
-\subfile{appendix/appendix_b}
+% \subfile{appendix/appendix_b}
-\subfile{appendix/appendix_c}
+% \subfile{appendix/appendix_c}
+% \subfile{appendix/appendix_d}
 
 \onecolumngrid
 

latex/references/references.bib

@@ -97,8 +97,32 @@
     urldate = {2023-10-22}
 }
 
-@misc{penning:fig,
+@online{penning:fig,
-    key = {something},
+    author = {Ludwig-Maximilians-Universität München},
+    title = {Penning traps},
     url = {https://www.med.physik.uni-muenchen.de/research/nuclear-science/nuclear-masses/mlltrap/layout/traps/index.html},
-    note = {Configuration figure a}
+    urldate = {2023-10-23},
+    note = {Configuration of a Penning trap, figure a}
+}
+
+@online{github:repo,
+    author = {Cory Alexander Balaton and Janita Ovidie Sandtrøen Willumsen},
+    url = {https://github.uio.no/FYS3150-G2-2023/Project-3},
+    urldate = {2023-10-24},
+    note = {FYS3150 Project 3 repo}
+}
+
+@online{scalasca,
+    author = {M. Geimer and F. Wolf and B.J.N. Wylie and E. Abraham and D. Becker and B. Mohr},
+    title = {Scalasca},
+    url = {https://www.scalasca.org/scalasca/about/about.html},
+    urldate = {2023-10-24},
+    note = {Tool to support performance optimization of parallel programs, measuring and analyzing runtime behavior.}
+}
+
+@online{scorep,
+    title = {Score-P: the Scalable Performance Measurement Infrastructure for Parallel Codes},
+    url = {https://perftools.pages.jsc.fz-juelich.de/cicd/scorep/tags/latest/html/},
+    urldate = {2023-10-24},
+    note = {Tool suite for profiling and event tracing.}
 }

latex/sections/abstract.tex

@@ -4,7 +4,7 @@
 \begin{document}
 
 \begin{abstract}
-    Add an abstract for the project?  
+    We have studied the behavior of singly-charged Calcium ions ($\text{Ca}^{+}$), inside an ideal Penning trap. With a numerical approach, we studied the equations of motion by implementing the forward Euler method \(FE\) and the 4th order Runge-Kutta \(RK4\). We found that RK4 approximates the solution with smaller relative error than the relative error of FE. In addition, we evaluated the methods by their rate of convergence. We found that RK4 has a higher convergence rate at approx. $4.0$, compared to FE at approx. $1.4$. For particles interacting we explored angular frequencies, and amplitudes, of the time-dependent potential applied to the particles. We found that angular frequency in the range $\omega_{V} \in (1.3, 1.4)$ MHz is effective in pushing out particles, even for amplitude $f = 0.1$.
 \end{abstract}
 
 \end{document}

latex/sections/conclusion.tex

@@ -2,11 +2,10 @@
 \graphicspath{{\subfix{../images/}}}
 
 \begin{document}
-
 \section{Conclusion}
+We studied the movement of particles confined by an ideal Penning trap, where we used iterative methods to simulate the particle behavior. We included the magnetic and electric field of the Penning trap, in addition to simulating the particles behavior when interaction with other each other. When we introduced the interaction, the movement in both radial direction and z-direction changed. From a circular path, to a more elliptical path, where the particles initial condition determine how it is affecting other particles path. 
 
+We also compared iterative methods, with the analytical solution, and found that the forward Euler \(FE\) method result in an approximation with a large relative error compared to the relative error of the 4th order Runge-kutta \(RK4\) method. In addition, we also found that RK4 has a higher convergence rate at approx. $4.0$, compared to FE at approx. $1.4$. Which suggest RK4 reach the solution faster than what FE, however, when we increase the number of time steps both methods result in similar relative error. When the number of calculations increase, and the number of time steps is sufficient, FE can be the better choise to conserve computational resources.
 
-
+When we explored the particles behavior at angular frequencies $\omega_{V} \in (0.2, 2.5)$ MHz, we found that particles are pushed out of the Penning trap when the amplitude of the applied time-dependent potential increase. The amplitude $f = 0.7$ result in particles being pushed out at most of the range of angular frequencies, whereas an amplitude $f = 0.1$ result in particles being pushed in a more narrow range. Since particles are being pushed out when the amplitude is low, there is likely a resonance frequency at around $1.4$ MHz.
-
-
 \end{document}

latex/sections/introduction.tex

@@ -4,10 +4,12 @@
 \begin{document}
 %
 \section{Introduction}
-We are surrounded by matter, which are made up of elementary particles. In the field of physics we want to understand the properties of these particles, measure their physical quantities, not to mention explain the origin of mass \cite{britannica:2023:matter}. 
+We are surrounded by matter, which are made up of elementary particles. In the field of physics we want to understand the properties of these particles, and measure their physical quantities. We want to describe the particles such that we can explain the origin of mass \cite{britannica:2023:matter}. 
 
-However, to study a particle, it is necessary isolate and contain it over time. The Penning trap is a device, able to confine charged particles for a period of time. This concept was evolved from F. M. Penning's implementation of magnetic fiels to a vaccum gauge, and J. R. Pierce's work with electron beams, and put into practice by Hans Dehmelt. In 1973 Dehmelt and his group of researchers were able to confine a particle and store it over several months \cite{vogel:2018:ch1}. 
+However, to study a particle, it is necessary isolate and contain it over time. The Penning trap is a device, able to confine charged particles for a period of time. This concept was evolved from F. M. Penning's implementation of magnetic fields to a vaccuum gauge, and J. R. Pierce's work with electron beams, and put into practice by Hans Dehmelt. In 1973 Dehmelt and his group of researchers were able to confine a particle and store it over several months \cite{vogel:2018:ch1}. 
 
-In practice, a Penning trap is not easy to obtain, and an experiment is both time consuming and expensive. A numerical approach, allow us to study the effects of the Penning trap on a charged particle, without the cost. We can use ordinary differential equations to model the particle's movement, confined within an Penning trap. Our focus will be on an ideal Penning trap, where an electrostatic field confines the particle in z-direction, and a magnetic field confines it in the radial direction. We will use numerical methods to model a single particle, and study the particle motion in radial direction. In addition, we will model a system of particles, and study their motion both with and without particle interaction.
+In practice, a Penning trap is not easy to obtain, and an experiment would be both time consuming and expensive. A numerical approach, allow us to study the effects of the Penning trap on a charged particle, without the equipment and material cost. With ordinary differential equations (ODE) we can model the movement of particles, confined within a Penning trap. 
+
+We will study an ideal Penning trap, where an electrostatic field confines the particle in z-direction, and a magnetic field confines it in the radial direction. We will use numerical methods to model a single particle, and study the particle motion in radial direction. In addition, we will model a system of particles, and study their motion both with and without particle interaction.
 %
 \end{document}

latex/sections/methods.tex

@@ -2,15 +2,48 @@
 \graphicspath{{\subfix{../images/}}}
 
 \begin{document}
-\section{Methods}
+\section{Methods}\label{sec:methods}
-% Problem 1
+\subsection{Theoretical background}\label{sec:theoretical_background}
-When we study the Penning traps effect on a particle with a charge $q$, we need to consider the forces acting on the particle. We can use Newton's second law \eqref{eq:newton_second} to determine the position of a particle as a function of time. In addition, we introduce the Lorentz force \eqref{eq:lorentz_force}, which describes the force on the particle. The position can be described by %
+In this experiment we study a particle with a charge $q$. To do so, we need to consider the forces acting on that particle. These forces are given by the Penning trap's electric field and magnetic field, an illustration can be found in fig. \ref{fig:penning_trap}. 
+\begin{figure}
+    \centering
+    \includegraphics[width=0.6\linewidth]{images/penning_trap.pdf}
+    \caption{The basic configuration of a hyperbolic Penning trap \cite{penning:fig}. Electric potential $V_{0}$ is applied to the electrodes in end caps (top and bottom), and in the ring (middle). The magnetic field is denoted as $\mathbf{B}$, the distance from the center to the end caps are denoted as $z_{0}$, and the radial distance as $r_{0}$.}
+    \label{fig:penning_trap}
+\end{figure}
+The electric field $\mathbf{E}$ in our model is related to the electric potential $V$ through 
+\begin{equation}\label{eq:electric_field_potential}
+    \mathbf{E} = - \nabla V,
+\end{equation} %
 %
+where we define the electric potential as
+\begin{equation}\label{eq:electric_potential}
+    V(x, y, z) = \frac{V_{0}}{2 d^{2}} (2z^{2} - x^{2} - y^{2}).
+\end{equation} %
+%
+The characteristic dimension $d = \sqrt{z_{0}^{2} + r_{0}^{2} / 2}$ determine the scale of the region between the electrodes. The magnetic field is homogeneous and defined as
+\begin{equation}\label{eq:magnetic_field}
+    \mathbf{B} = B_{0} \hat{e}_{z} = (0, 0, B),
+\end{equation} %
+%
+where $B_{0} > 0$ determines the strength of the field. The electric potential $V_{0}$ is applied to the electrodes, where the end caps are positively charged and the ring is negatively charged. The particle is confined by the electric field in the z-direction, however, it is not confined in the radial direction (xy-plane). The magnetic field is necessary to ensure the particle is fully confined in the Penning trap, and will force the particle to move in a circular orbit. 
+
+First, we consider Newton's second law \eqref{eq:newton_second}, to determine the position of the particle. 
+\begin{equation}\label{eq:newton_second}
+    m \ddot{\mathbf{r}} = \sum_{i} \mathbf{F}_{i}
+\end{equation} % 
+%
+In addition, we introduce the Lorentz force \eqref{eq:lorentz_force}, which describes the force acting on the particle.
+\begin{equation}\label{eq:lorentz_force}
+    \mathbf{F} = q \mathbf{E} + q \mathbf{v} \times \mathbf{B},
+\end{equation} %
+%
+where $\ddot{\mathbf{r}} = d^{2} \mathbf{r} / dt^{2}$. We combine eq. \eqref{eq:newton_second} and eq. \eqref{eq:lorentz_force} in
 \begin{equation}\label{eq:newton_lorentz}
-    m \ddot{\mathbf{r}} = (q \mathbf{E} + q \mathbf{v} \times \mathbf{B}).
+    m \ddot{\mathbf{r}} = (q \mathbf{E} + q \mathbf{v} \times \mathbf{B}),
-\end{equation}
+\end{equation} %
 %
-Using eq. \eqref{eq:newton_lorentz} we derive the differential equations in sec. \ref{sec:eq_motion}, and can rewrite them as
+and derive the differential equations in appendix \ref{sec:eq_motion}. We can rewrite and simplify these equations as
 \begin{align}
     \label{eq:motion_x}
     \ddot{x} - \omega_{0} \dot{y} - \frac{1}{2} \omega_{z}^{2} x &= 0, \\
@@ -20,32 +53,36 @@ Using eq. \eqref{eq:newton_lorentz} we derive the differential equations in sec.
     \ddot{z} + \omega_{z}^{2} z &= 0,
 \end{align} %
 %
+where
+\begin{equation*}
+    \omega_{0} \equiv \frac{q B_{0}}{m}, \quad \omega_{z}^{2} \equiv \frac{2 q V_{0}}{m d^{2}}.
+\end{equation*} %
+%
 We find the general solution for eq. \eqref{eq:motion_z}
 \begin{equation}\label{eq:eq_general}
     z(t) = c_{1} e^{i \omega_{z} t} + c_{2} e^{-i \omega_{z} t},
 \end{equation} %
-derived in sec. \ref{sec:eq_general}. Continuing, we will use a Calcium ion with a single positive charge. That is, we assume the charge of the particle is $q > 0$.
+derived in appendix \ref{sec:eq_general}. Continuing, we will use a Calcium ion with a single positive charge. That is, we assume the charge of the particle is $q > 0$.
 
-% Problem 2
+Since eq. \eqref{eq:motion_x} and eq. \eqref{eq:motion_y} are coupled, we want to rewrite them as a single differential equation. We derive this in appendix \ref{sec:eq_complex}, and the resulting equation is given by
-Since eq. \eqref{eq:motion_x} and eq. \eqref{eq:motion_y} are coupled, we want to rewrite them as a single differential equation. We derive this in sec. \ref{sec:eq_complex}, and the resulting equation is given by
 \begin{equation}\label{eq:single_differential}
     \ddot{f} + i \omega_{0} \dot{f} - \frac{1}{2} \omega_{z}^{2} f = 0. 
-\end{equation}
+\end{equation} %
-
+%
-% Problem 3
 Eq. \eqref{eq:single_differential} has a general solution given by
 \begin{equation}\label{eq:general_solution}
     f(t) = A_{+}e^{-i(\omega_{+} t + \phi_{+})} + A_{-}e^{-i(\omega_{-} t + \phi_{-})}.
 \end{equation} %
-The amplitude $A_{+}$ and $A_{-}$ are positive, the phases $\phi_{+}$ and $\phi_{-}$ are constant, and the rate is given by
+%
-\begin{equation*}
+The amplitude $A_{+}$ and $A_{-}$ are positive, the phases $\phi_{+}$ and $\phi_{-}$ are constant, and the angular rate is given by
+\begin{equation}\label{eq:angular_rate}
     \omega_{\pm} = \frac{\omega_{0} \pm \sqrt{\omega_{0}^{2} - 2 \omega_{z}^{2}}}{2}.
-\end{equation*}
+\end{equation} 
 We find the physical coordinates at a given time $t$ using 
 \begin{equation*}
     x(t) = \text{Re}f(t), \quad y(t) = \text{Im}f(t),
 \end{equation*} %
-and eq. \eqref{eq:general_solution}. We can rearrange the right hand side of the derived expression in \ref{sec:eq_coord}, and find the physical coordinates
+and eq. \eqref{eq:general_solution}. We can rearrange the right hand side of the expression derived in \ref{sec:eq_coord}, to find the physical coordinates
 \begin{align}\label{eq:physical_coord}
     x(t) &= A_{+} \cos(\omega_{+} t + \phi_{+}) + A_{-} \cos(\omega_{-} t + \phi_{-}) \\
     y(t) &= - A_{+} \sin(\omega_{+} t + \phi_{+}) - A_{-} \sin(\omega_{-} t + \phi_{-})
@@ -55,104 +92,74 @@ However, to obtain a bound on the particle's movement in the radial direction, w
 \begin{align*}
     |f(t)| &= \sqrt{(x(t))^{2} + (y(t))^{2}}. \\
 \end{align*} %
+Which means we have to put constraint on the values of $\omega_{0}$ and $\omega_{z}^{2}$ in order to find a bounded solution $|f(t)| < \infty$. We derive an expression for this in appendix \ref{sec:bounded_solution}. Since we have assumed the particle charge $q > 0$, and the mass $m > 0$, the constraints are
+\begin{equation}
+    \frac{q}{m} > \frac{4 V_{0}}{d^{2} B_{0}^{2}}.
+\end{equation} %
+% 
 When $t \rightarrow \infty$, the upper and lower limits are
 \begin{align}
-    R_{+} &= A_{+} + A_{-} \\
+    \label{eq:upper_b}
+    R_{+} &= A_{+} + A_{-}, \\
+    \label{eq:lower_b}
     R_{-} &= |A_{+} - A_{-}|, 
 \end{align}
-derived in sec. \ref{sec:bounded_solution}.
+which we derive in appendix \ref{sec:upper_lower_bound}.
 
-% Problem 3
+In this experiment, we will use the initial conditions for the particle given in table \ref{tab:initial_particle_cond}. From eq. \eqref{eq:motion_z} and the initial conditions, we find the specific solution
-To obtain the bounded solution, we need to consider the expression
 \begin{equation*}
-    \alpha = (\omega_{+} - \omega_{-}) t +( \phi_{+} -  \phi_{-}).
+    z(t) = z_{0} \cos (\omega_{z} t)
 \end{equation*}
-When $t \rightarrow \infty$ the constant phases will not affect the expression, we obtain a bounded solution if we consider the rate such that $|f(t)| < \infty$. That is, we need $\omega_{0}^{2} - 2 \omega_{z}^{2}$ in eq. \eqref{eq:angular_rate} to be real.
+For the particle movement in the radial direction, we find the specific solution using eq. \eqref{eq:general_solution} with 
 \begin{align*}
-    \omega_{0}^{2} - 2 \omega_{z}^{2} = 0, \quad \omega_{0} > 2 \omega_{z}
+    A_{+} &= \frac{v_{0} + \omega_{-} x{0}}{\omega_{-} - \omega_{+}}, & A_{-} &= - \frac{v_{0} + \omega_{+} x{0}}{\omega_{-} - \omega_{+}}, \\
+    \phi_{+} &= 0, & \phi_{-} &= 0.
 \end{align*}
-\begin{equation}\label{eq:angular_rate}
+
-    \omega_{\pm} = \frac{\omega_{0} \pm \sqrt{\omega_{0}^{2} - 2 \omega_{z}^{2}}}{2}
+We will also consider the numerical simulation of multiple particles, confined in a Penning trap. The particles will experience a repelling force from each other, giving the electric field 
+\begin{equation}\label{eq:electric_field_interaction}
+    \mathbf{E} = k_{e} \sum_{j = 1}^{n} q_j \frac{\mathbf{r} - \mathbf{r}_{j}}{|\mathbf{r}_{i} - \mathbf{r}_{j}|^{3}}.
+\end{equation} %
+The Coulomb constant $k_{e}$ value can be found in table \ref{tab:constants}. 
+
+For multiple particles we have to modify the equations of motion, by adding a term for the force a given particle experience at a given point. When we scale eq. \eqref{eq:electric_field_interaction} by charge and mass, we get a new set of equations of motion
+\begin{align}
+    \label{eq:coulomb_motion_x}
+    \ddot{x}_{i} - \omega_{0,i} \dot{y}_{i} - \frac{1}{2} \omega_{z,i}^{2} x_{i} - k_{e} \frac{q_{i}}{m_{i}}\sum_{j \neq i} q_j \frac{x_{i} - x_{j}}{|\mathbf{r}_{i} - \mathbf{r}_{j}|^{3}} &= 0, \\
+    \label{eq:coulomb_motion_y}
+    \ddot{y}_{i} - \omega_{0,i} \dot{x}_{i} - \frac{1}{2} \omega_{z,i}^{2} y_{i} - k_{e} \frac{q_{i}}{m_{i}}\sum_{j \neq i} q_j \frac{y_{i} - y_{j}}{|\mathbf{r}_{i} - \mathbf{r}_{j}|^{3}} &= 0, \\
+    \label{eq:coulomb_motion_z}
+    \ddot{z}_{i} + \omega_{z,i}^{2} z_{i} - k_{e} \frac{q_{i}}{m_{i}}\sum_{j \neq i} q_j \frac{z_{i} - z_{j}}{|\mathbf{r}_{i} - \mathbf{r}_{j}|^{3}} &= 0,
+\end{align} %
+where $i$ and $j$ denote the particle indices. When we include a time-dependence to the applied potential we make a replacement of the initial electric potential
+\begin{equation}\label{eq:pertubation}
+    V_{0} \rightarrow V_{0} (1 + f \cos (\omega_{V} t)),
 \end{equation}
+where $f$ denotes the amplitude and $\omega_{V}$ the angular rate.
 
-We find the solution 
+\subsection{Algorithms and implementation}\label{sec:algo_implementation}
-Physical properties given by newtons second law \eqref{eq:newton_second}
+When we consider a multi-particle system, we adapt the notation in our derived equations. The force acting on a particle in the Penning trap, which experience particle interaction, can be written as
-
+\begin{equation}\label{eq:lorentz_force_interaction}
-The particle moves and its position can be determined using newton. where the electric field
+    \mathbf{F}_{i}(t, \mathbf{v}_{i}, \mathbf{r}_{i}) = q_{i} \mathbf{E}(t, \mathbf{r}_{i}) + q_{i} \mathbf{v}_{i} \cross \mathbf{B} - \mathbf{E}_{p}(t, \mathbf{r}_{i}),
-
-
-
-
-
-\begin{table}[H]
-    \centering
-        \begin{tabular}[c]{lll}
-            Constant & Value &  \\
-            \hline
-            $B_{0}$ & $1.00 T$ & $9.65 \cross 10^{1} \frac{u}{(\mu s) e}$ \\
-            $V_{0}$ & $25.0 mV$ & $2.41 \cross 10^{6} \frac{u (\mu m)^{2}}{(\mu s)^{2} e}$ \\
-            $d$ & $500 \ \mu m$ & \\
-            \hline
-        \end{tabular}
-    \caption{Default configuration of the Penning trap, where the value of T and V can be found in table \ref{tab:constants}.}
-    \label{tab:penning_config}
-\end{table}
-
-\subsection{Dealing with a multi--particle system}
-
-For a multi-particles system, we need to modify $\vb{F}$ to account for the 
-force of other particles in the system acting upon each other. To do that, we
-add another term to $\vb{F}$
-\begin{equation}
-    \vb{F}_i(t, \vb{v}_i, \vb{r}_i) = q_i \vb{E}(t, \vb{r}_i) + q_i \vb{v}_i \cross \vb{B} - \vb{E}_p(t, \vb{r}_i),
-\end{equation}
-where $i$ and $j$ are particle indices and
-\begin{equation}
-    \vb{E}_p(t, \vb{r}_i) = q_i k_e \sum_{j \neq i} 
-    q_j \frac{\vb{r_i} - \vb{r_j}}{\left| \vb{r_i} - \vb{r_j} \right|^3}.
-    \label{eq:}
 \end{equation}
+where $i$ and $j$ still denotes the particle indices, and $\mathbf{E}_{p}$ denotes the force exerted on particle $i$ by particle $j$.
+\begin{equation}\label{eq:electric_force_interaction}
+    \mathbf{E}_{p}(t, \mathbf{r}_{i}) = k_{e} q_{i} \sum_{j \neq i} 
+    q_{j} \frac{\mathbf{r}_{i} - \mathbf{r}_{j}}{| \mathbf{r}_{i} - \mathbf{r}_{j} |^{3}}.
+\end{equation} %
+In addition, we adapt the notation of \eqref{eq:newton_second}, and define the first derivative of the particle position as
+\begin{align}
+    \label{eq:r_derivative}
+    \frac{d \mathbf{r}_{i}}{dt} = \dot{\mathbf{r}}_{i} &= \mathbf{v}_{i}, \\
+    \label{eq:v_derivative}
+    \frac{d \mathbf{v}_{i}}{dt} = \dot{\mathbf{v}}_{i} &= \frac{\mathbf{F}_{i} (t, \mathbf{v}_{i}, \mathbf{r}_{i}}{m_{i}},
+\end{align}
 
-Newton's second law for a particle $i$ is then
+We first implemented the forward Euler method, using the expression for a coupled system given in appendix \ref{sec:algo_euler}. We define the forward Euler algorithm in \ref{algo:forward_euler}. 
-\begin{equation}
-    \frac{d^2\vb{r}_i}{dt^2} = \frac{\vb{F}_i\left(t, \frac{d\vb{r}_i}{dt}, \vb{r_i}\right)}{m_i},
-    \label{eq:newtonlaw2}
-\end{equation}
-We can then rewrite the second order ODE from equation~\ref{eq:newtonlaw2} 
-into a set of coupled first order ODEs. 
-We now rewrite Newton's second law of motion as
-\begin{equation}
-    \begin{split}
-        \frac{d\vb{r}_i}{dt} &= \vb{v}_i \\
-        \frac{d\vb{v}_i}{dt} &= \frac{\vb{F}_i(t, \vb{v}_i, \vb{r}_i)}{m_i}.
-    \end{split}
-    \label{eq:coupled}
-\end{equation}
-
-\subsection{Forward Euler}
-
-For a particle $i$, the forward Euler method for a coupled system is 
-expressed as
-\begin{equation}
-    \begin{split}
-        \vb{r}_{i,j+1} &= \vb{r}_{i,j} + h \cdot \frac{d\vb{r}_{i,j}}{dt} = \vb{r}_{i,j} + h \cdot \vb{v}_{i,j} \\
-        \vb{v}_{i,j+1} &= \vb{v}_{i,j} + h \cdot \frac{\vb{v}_{i,j}}{dt} = \vb{v}_{i,j} + h \cdot \frac{\vb{F}\left( t_{j}, \vb{v}_{i,j}, \vb{r}_{i,j} \right)}{m},
-    \end{split}
-\end{equation}
-for particle $i$ where $j$ is the current time step of the particle, 
-$m$ is the mass of the particle, and $h$ is the step length.
-
-When dealing with a multi-particle system, we need to ensure that we do not
-update the position of any particles until every particle has calculated their
-next step. An easy way of doing this is to create a copy of all the particles,
-then update the copy, and when all the particles have calculated their next 
-step, simply replace the particles with the copies. 
-Algorithm~\ref{algo:forwardeuler} provides an overview on how that can be achieved. % Make this better
-
 \begin{figure}[H]
     \begin{algorithm}[H]
     \caption{Forward Euler method}
-    \label{algo:forwardeuler}
+    \label{algo:forward_euler}
         \begin{algorithmic}
             \Procedure{Evolve forward Euler}{$particles, dt$}
             \State $N \leftarrow \text{Number of particles in } particles$
@@ -166,56 +173,7 @@ Algorithm~\ref{algo:forwardeuler} provides an overview on how that can be achiev
     \end{algorithm}
 \end{figure}
 
-\subsection{4th order Runge-Kutta}
+We also implemented the 4th order Runge-Kutta (RK4) method, using the expression given in appendix \ref{sec:algo_rk4}. We define the RK4 algorithm in \ref{algo:rk4}. $\mathbf{F}$ does not take any arguments, however, the total force acting on the particle is calculated using the value of position and velocity within $particles$. 
-
-For a particle $i$, we can express the 4th order Runge-Kutta (RK4) method as
-\begin{equation}
-    \begin{split}
-        \vb{v}_{i,j+1} &= \vb{v}_{i,j} + \frac{h}{6} \left( \vb{k}_{\vb{v},1,i} 
-            + 2\vb{k}_{\vb{v},2,i} + 2\vb{k}_{\vb{v},3,i} + \vb{k}_{\vb{v},4,i} 
-            \right) \\
-        \vb{r}_{i,j+1} &= \vb{r}_{i,j} + \frac{h}{6} \left( \vb{k}_{\vb{r},1,i} 
-            + 2\vb{k}_{\vb{r},2,i} + 2\vb{k}_{\vb{r},3,i} + \vb{k}_{\vb{r},4,i} 
-            \right),
-    \end{split}
-\end{equation}
-where
-\begin{equation}
-    \begin{split}
-        \vb{k}_{\vb{v},1,i} &= \frac{\vb{F}_i(t_j, \vb{v}_{i,j}, 
-            \vb{r}_{i,j})}{m} \\
-        \vb{k}_{\vb{r},1,i} &= \vb{v}_{i,j} \\
-        \vb{k}_{\vb{v},2,i} &= \frac{\vb{F}_i(t_j+\frac{h}{2}, \vb{v}_{i,j} 
-            + h \cdot \frac{\vb{k}_{\vb{v},1,i}}{2}, \vb{r}_{i,j} 
-            + h \cdot \frac{\vb{k}_{\vb{r},1,i}}{2})}{m} \\
-        \vb{k}_{\vb{r},2,i} &= \vb{v}_{i,j} 
-            + h \cdot \frac{\vb{k}_{\vb{v},1,i}}{2} \\
-        \vb{k}_{\vb{v},3,i} &= \frac{\vb{F}_i(t_j+\frac{h}{2}, \vb{v}_{i,j} 
-            + h \cdot \frac{\vb{k}_{\vb{v},2,i}}{2}, \vb{r}_{i,j} 
-            + h \frac{\cdot \vb{k}_{\vb{r},2,i}}{2})}{m} \\
-        \vb{k}_{\vb{r},3,i} &= \vb{v}_{i,j} 
-            + h \cdot \frac{\vb{k}_{\vb{v},2,i}}{2} \\
-        \vb{k}_{\vb{v},4,i} &= \frac{\vb{F}_i(t_j+h, \vb{v}_{i,j} 
-            + h \cdot \vb{k}_{\vb{v},3,i}, \vb{r}_{i,j} 
-            + h \cdot \vb{k}_{\vb{r},3,i})}{m} \\
-        \vb{k}_{\vb{r},4,i} &= \vb{v}_{i,j} 
-            + h \cdot \frac{\vb{k}_{\vb{v},1,i}}{2}.
-    \end{split}
-\end{equation}
-
-In order to find each $\vb{k}_{\vb{r},i}$ and $\vb{k}_{\vb{v},i}$, 
-we need to first compute all $\vb{k}_{\vb{r},i}$ and $\vb{k}_{\vb{v},i}$ 
-for all particles, then we can update the particles in order to compute 
-$\vb{k}_{\vb{r},i+1}$ and $\vb{k}_{\vb{v},i+1}$. In order for the algorithm
-to work, we need to save a copy of each particle before starting so that we
-can update the particles correctly for each step.
-
-This approach would require 8 loops to be able to complete the calculation since
-we cannot update the particles until after all $\vb{k}$ values have been 
-computed, however if we create a temporary array that holds particles, we can
-put the updated particles in there, and then use that array in the next loop,
-and would reduce the required amount of loops down to 4.
-
 \begin{figure}
     \begin{algorithm}[H]
     \caption{RK4 method}
@@ -267,12 +225,12 @@ and would reduce the required amount of loops down to 4.
             \State $\vb{k}_{\vb{v},4,i} \leftarrow \frac{\vb{F}}{m}$
 
             \State $tmp\_p_i.\vb{r} \leftarrow orig\_p_i.\vb{r} + \frac{dt}{6} 
-                \cdot \left( \vb{k}_{\vb{r},1,i} + \vb{k}_{\vb{r},2,i} 
+                \cdot \left( \vb{k}_{\vb{r},1,i} + 2 \vb{k}_{\vb{r},2,i} 
-                + \vb{k}_{\vb{r},3,i} + \vb{k}_{\vb{r},4,i} \right)$
+                + 2 \vb{k}_{\vb{r},3,i} + \vb{k}_{\vb{r},4,i} \right)$
 
             \State $tmp\_p_i.\vb{v} \leftarrow orig\_p_i.\vb{v} + \frac{dt}{6} 
-                \cdot \left( \vb{k}_{\vb{v},1,i} + \vb{k}_{\vb{v},2,i} 
+                \cdot \left( \vb{k}_{\vb{v},1,i} + 2 \vb{k}_{\vb{v},2,i} 
-                + \vb{k}_{\vb{v},3,i} + \vb{k}_{\vb{v},4,i} \right)$
+                + 2 \vb{k}_{\vb{v},3,i} + \vb{k}_{\vb{v},4,i} \right)$
             
             \EndFor
 
@@ -280,17 +238,20 @@ and would reduce the required amount of loops down to 4.
             \EndProcedure
         \end{algorithmic}
     \end{algorithm}
-    $\vb{F}$ in the algorithm does not take any arguments as it uses the 
-    velocities and positions of the particles inside the array $particles$ to
-    calculate the total force acting on particle $i$.
 \end{figure}
+For a system of particles, the positional value of one particle affects the rest of the particles. We calculated the next step for all particles, updated their values in a copy, and updated the positional value based on the original values. We simulated the movement of particles confined in a Penning trap. All simulations used the initial conditions for particle 1 and 2 given in table \ref{tab:initial_particle_cond}. 
  
-
+\subsection{Testing and error analysis}
-\subsection{Testing the simulation}
+We implemented a test suite to validate our results, details can be found in appendix \ref{sec:testing}. To analyze our results, we compared the relative error of the analytical solution with the result from the implemented forward Euler method and RK4 method. We simulated four times using time steps given in table \ref{tab:time_steps}, and estimated the error convergence rate $r_{\text{err}}$ for both methods using 
-
+\begin{equation}
-\subsection{Relative error and error convergance rate}
+    r_{\text{err}} = \frac{1}{3} \sum_{k=2}^{4} \frac{\log(\Delta_{\text{max},k} / \Delta_{\text{max}, k-1})}{\log(h_{k} / h_{k-1})}. 
+\end{equation}
+The maximum error of a simulation $k$ with a step size $h_{k}$ is given by
+\begin{equation*}
+    \Delta_{\text{max},k} = \max_{i} |\mathbf{r}_{i, \text{exact}} - \mathbf{r}_{i}|,
+\end{equation*}
+where $\mathbf{r}_{i, \text{exact}}$ is the analytical solution, and $\mathbf{r}_{i}$ is the computed result.
 
 \subsection{Tools}
-The numerical methods implemented in C++, are parallelized using OpenMP \cite{openmp:2018}. We used the Python library matplotlib \cite{hunter:2007:matplotlib} to produce all the plots, and seaborn \cite{waskom:2021:seaborn} to set the theme in the figures.
+The numerical methods are implemented in C++, and parallelized using \verb|OpenMP| \cite{openmp:2018}. In addition, we used the profiler \verb|Scalasca| \cite{scalasca} and \verb|Score-P| \cite{scorep} event tracking, during development. We used 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}

latex/sections/results.tex

@@ -2,88 +2,122 @@
 \graphicspath{{\subfix{../images/}}}
 
 \begin{document}
-
 \section{Results and Discussion}
-% Single particle
+The simulations were performed using a constant configuration for the Penning trap, where the values of $B_{0}$, $V_{0}$, and $d$ can be found in table \ref{tab:penning_config}. We also used a constant configuration for the particles, found in table \ref{tab:particle_config}. Initial position of particle 1 $p_{1}$ was set to $(20, 0, 20) \mu m$ with velocity $(0, 25, 0) \mu m/ \mu s$, whereas the position of particle 2 $p_{2}$ was set to $(25, 25, 0) \mu m$ with velocity $(0, 40, 5) \mu m/ \mu s$.
-We simulated the movement of particles confined in a Penning trap. All simulations used the initial conditions for particle 1 and 2 given in table \ref{tab:initial_condition_particles}. 
 
-First we simulated a single particle for $50 \mu s$, approximating the particle's motion using the RK4 method. In addition we compared the motion of particle 1 with the analytical solution in figure \ref{fig:single_particle}. What we see is a complete overlap of the analytical solution completely overlap the approximated, suggest that the simulation result is good.
+First, we simulated a single particle for $50 \mu s$, approximating the particle's motion using the RK4 method. In addition we compared the result with the analytical solution in figure \ref{fig:single_particle}. The approximated solution completely overlap the analytical, suggesting the implemented method is approximating the solution with minimal error. The angular frequency is
-% Add something about why the simulated result is good, cos(wt) when w is
+\begin{align*}
-% 
+    \omega_{z} &= \sqrt{\frac{2 q V_{0}}{m d^{2}}} = \sqrt{\frac{2 \cdot 1 \cdot 9.65}{40.078}} \approx 0.694 \text{ rad/\textmu s},
-\begin{table}[H]
+\end{align*}
-    \centering
+which result in a period (P) of
-        \begin{tabular}[c]{lll}
+\begin{align*}
-            Particle & Position & Velocity \\
+    \text{P} &= \frac{2 \pi}{|\omega_{z}|} \approx 9.054 \text{\textmu s}.
-            \hline
+\end{align*}
-            $p_{1}$ & $(20, 0, 20) \mu m$ & $(0, 25, 0) \mu m/ \mu s$ \\
+From figure \ref{fig:single_particle} we see that the period, the time it takes for the particle to reach the same position, is close to 9 \textmu s. This is what we would expect given the value of $\omega_{z}$.
-            $p_{2}$ & $(25, 25, 0) \mu m$ & $(0, 40, 5) \mu m/ \mu s$ \\
+% Figure: single_particle.pdf
-            \hline
-        \end{tabular}
-    \caption{Initial position and velocity of particle 1 ($p_{1}$) particle 2 ($p_{2}$), where the analytical solution is given by $z(t) = z_{0} \cos (\omega_{z} t)$}
-    \label{tab:initial_condition_particles}
-\end{table}
 \begin{figure}[H]
     \centering
     \includegraphics[width=\linewidth]{images/single_particle.pdf}
-    \caption{A single particle in the Penning trap, approximated and analytical motion in z-direction.}
+    \caption{Movement of a single particle in z-direction. Approximated using the 4th order Runge-Kutta method, compared to the analytical solution $z(t) = z_{0} \cos (\omega_{z} t)$.}
     \label{fig:single_particle}
 \end{figure}
-% Add equations of motion for particle with interaction eq. (18, 19, 20)
+
-% Multiple particles
+To evaluate the implemented methods we simulated the particle using both forward Euler (FE) and RK4, with different time steps given in \ref{tab:time_steps}. Again, we simulated the particle movement for $50$ \textmu s, and estimated the relative error of each method. The result can be found in figure \ref{fig:relative_error}. The relative error of FE is large compared to the relative error of RK4. The error convergence rate for FE is $r_{\text{err}} \approx 1.39652$, and RK4 $r_{\text{err}} \approx 3.99998$, which means RK4 converges faster compared to FE.
-% Add initial condition of Penning trap
+\begin{figure}[H]
-We will now consider the Penning trap with initial conditions given in table \ref{tab:initial_condition_penning}, and simulate using one or two particles. In addition, we simulate two particles both with and without interactions, the result is found in figure \ref{fig:two_particles}. When we add interaction between the particles, they both still follow the same inherent path. However, we observe a small shift in both particle's movement.
-\begin{table}[H]
     \centering
-    \begin{tabular}{lll}
+    \includegraphics[width=\linewidth]{images/relative_error.pdf}
-        $B_0$ & $V_{0}$ & $d$ \\
+    \caption{Relative error measured for 4th order Runge-Kutta \(a\) and foreward Euler \(b\). Important to notice when reading the plot, that the scale of measured relative error differ between the methods.}
-        \hline
+    \label{fig:relative_error}
-    \end{tabular}
-    \caption{Caption}
-    \label{tab:initial_condition_penning}
-\end{table}
-%
-\begin{figure}
-    \centering
-    \includegraphics[width=\linewidth]{images/plot_2_particles_xy.pdf}
-    \caption{Movement of two particles in the xy-plane. $\hat{p}_{1}$ and $\hat{p}_{2}$ include particle interaction, whereas $p_{1}$ and $p_{2}$ does not include particle interaction.}
-    \label{fig:two_particles}
 \end{figure}
-% Phase space plot
+RK4 is a more complex method, which requires four times the calculations compared to FE. When the simulation is done with $n_{4} = 32000$ time steps, the difference in relative error decreases. The cost of calculation exceed the benefit of accuracy. That is, at a large number of time steps we could find a satisfying approximation to the solution using FE, with less calculations. 
-When we simulate two particles, we can see the effect of interaction between the particles in the xy-plane in fig. \ref{fig:phase_space_2x} and in the z-direction in fig. \ref{fig:phase_space_2z}. What we observe is a very small shift in position for particle 1 in x-direction, whereas particle 2 does not have a visible shift. In the z-direction, however, the oscillation of particle 2 experience a greater shift. Particle 2 experience the force of particle 1 such that particle 2 moves larger distance.
+
-%
+Next, we simulated the movement of two particles in radial direction. In figure \ref{fig:two_particles_radial}, we see the path of both particles without interaction. The starting position of particle 1 $p_{1}$ is closer to the center of the Penning trap, than the starting position of $p_{}$. Both particles move in an orbital path arount the center, however, the starting velocity of each particles determine the minimum and maximum distance a particle moves from the center. We find these values bounds using eq. \ref{eq:upper_b} and eq. \ref{eq:lower_b}
-\begin{figure}
+\begin{align*}
+    R_{+} &= A_{+} + A_{-} \\
+    &= -12.3232 \text{\textmu m} + 32.3232 \text{\textmu m} \\
+    &= 20.00 \text{\textmu m},
+\end{align*}
+and
+\begin{align*}
+    R_{-} &= |A_{+} - A_{-}| \\
+    &= |-12.3232 \text{\textmu m} - 32.3232 \text{\textmu m}| \\
+    &= 44.6464 \text{\textmu m}.
+\end{align*}
+The distance is similar that in figure \ref{fig:two_particles_radial}. In addition, we can see from \ref{fig:phase_nointer_x} that the distance differ between $p_{1}$ and $p_{2}$. When we study the particle movement in z-direction, in figure \ref{fig:phase_nointer_z}, we see a circular movement. Again, the distance the particle moves from the center is determined by its initial conditions. 
+
+In figure \ref{fig:two_particles_radial_interaction} we see the movement in radial direction, where particle interaction is included. The orbital path is similar to what we see without particle interaction. However, the additional force acting on the particle is affecting it's trajectory. The distance each particle is moving from the center, differs depending on the position. Which is also seen in figure \ref{fig:phase_inter_x}, where we see the trajectory of the particle at a given position. When we study the phase space plot, where particle interaction is included, in figure \ref{fig:phase_inter_z} a similar difference is observed. In figure \ref{fig:3d_particles} we see the movement of two particles in both the radial direction and z-direction. The figure include particle movement with and without interaction, where we can observe the change in movement.% Something about the trajectory related to the interaction?
+
+\begin{figure}[H]
     \centering
-    \includegraphics[width=\linewidth]{images/phase_space_no_interaction_x.pdf}
+    \includegraphics[width=\linewidth]{images/simulate_2_particles_no_interaction_xy.pdf}
-    \caption{Phase space plot of two particles in x-direction.}
+    \caption{Path of two particles, $p_{1}$ and $p_{2}$, without interaction. The black dot marks the starting position of each particle.}
-    \label{fig:phase_space_2x}
+    \label{fig:two_particles_radial}
-\end{figure}
-%
-\begin{figure}
-    \centering
-    \includegraphics[width=\linewidth]{images/phase_space_interaction_z.pdf}
-    \caption{Phase space plot of two particles in z-direction.}
-    \label{fig:phase_space_2z}
 \end{figure}
 
-\begin{figure}
+\begin{figure}[H]
+    \centering
+    \includegraphics[width=\linewidth]{images/simulate_2_particles_interaction_xy.pdf}
+    \caption{Path of two particles, $p_{1}$ and $p_{2}$, with interaction. The black dot marks the starting position of each particle.}
+    \label{fig:two_particles_radial_interaction}
+\end{figure}
+
+\begin{figure}[H]
+    \centering
+    \includegraphics[width=\linewidth]{images/phase_space_no_interaction_x.pdf}
+    \caption{Particle trajectory at a given position $x$, without particle interaction. The black dot marks the starting position of each particle.}
+    \label{fig:phase_nointer_x}
+\end{figure}
+
+\begin{figure}[H]
+    \centering
+    \includegraphics[width=\linewidth]{images/phase_space_no_interaction_z.pdf}
+    \caption{Particle trajectory at a given position $z$, without particle interaction. The black dot marks the starting position of each particle.}
+    \label{fig:phase_nointer_z}
+\end{figure}
+
+\begin{figure}[H]
+    \centering
+    \includegraphics[width=\linewidth]{images/phase_space_interaction_x.pdf}
+    \caption{Particle trajectory at a given position $x$, with particle interaction. The black dot marks the starting position of each particle.}
+    \label{fig:phase_inter_x}
+\end{figure}
+
+\begin{figure}[H]
+    \centering
+    \includegraphics[width=\linewidth]{images/phase_space_interaction_z.pdf}
+    \caption{Particle trajectory at a given position $z$, with particle interaction. The black dot marks the starting position of each particle.}
+    \label{fig:phase_inter_z}
+\end{figure}
+
+\begin{figure}[H]
     \centering
     \includegraphics[width=\linewidth]{images/3d_plot.pdf}
-    \caption{3D plot of particles-}
+    \caption{The movement of particles in a Penning trap. Where $p_{1}$ and $p_{2}$ is without interaction, $\hat{p}_{1}$ and $\hat{p}_{2}$ is without interaction.}
     \label{fig:3d_particles}
 \end{figure}
 
-\begin{figure}
+Finally, by subjecting the system to a time-dependent field, making the replacement in \ref{eq:pertubation}, we study the fraction of particles left at different amplitudes $f$. We can see how the different amplitudes lead to loss of particles, at different angular frequencies $\omega_{V}$ in \ref{fig:wide_sweep}. We study frequencies in the range $\omega_{V} \in (0.2, 2.5)$ MHz, with steps of $0.02$ MHz, and find that angular frequencies in the range $(1.0, 2.5)$ is effective in pushing the particles out of the Penning trap. 
-    \centering
-    \includegraphics[width=\linewidth]{images/relative_error.pdf}
-    \caption{Relative error of RK4 and forward Euler method.}
-    \label{fig:rel_err}
-\end{figure}
 
+We explore the range $\omega_{V} \in (1.0, 1.7)$ MHz closer in figure \ref{fig:narrow_sweep}, and observe a loss of particles for amplitude $f_{1} = 0.1$ in the range $\omega_{V} \in (1.3, 1.5)$, for $f_{2} = 0.4$ in the range $\omega_{V} \in (1.2, 1.6)$, and $f_{3} = 0.7$ in the entire range explored. When we introduce particle interaction, we observe the amplitude with higher value will result in a larger bound for the particle movement, and particles are easily pushed out. When we add particle interaction the angular frequencies, which result in particles being pushed out, are similar as in no interaction. As we see in figure \ref{fig:narrow_sweep_interactions} the particle's behavior, when interactions are added, is disrupted and add to the force resulting in the particle being pushed out of the Penning trap. We see that fewer particles tend to be left in the Penning trap, than without any particle interactions. We also see that the fraction of particles left is more unpredictable with particle interactions. 
 
-\begin{figure}
+\begin{figure}[H]
     \centering
     \includegraphics[width=\linewidth]{images/particles_left_wide_sweep.pdf}
-    \caption{Fraction of particles left in the Penning trap, with a given amplitude $f$.}
+    \caption{Exploring the behavior of particles, where the amplitude of time-dependent potential $f = [0.1, 0.4, 0.7]$, for angular frequency $\omega_{V} \in (0.2, 2.5)$ MHz.}
-    \label{fig:particles_left}
+    \label{fig:wide_sweep}
+\end{figure}
+
+\begin{figure}[H]
+    \centering
+    \includegraphics[width=\linewidth]{images/particles_left_narrow_sweep.pdf}
+    \caption{Exploring the behavior of particles, where the amplitude of time-dependent potential $f = [0.1, 0.4, 0.7]$, for angular frequency $\omega_{V} \in (1.1, 1.7)$ MHz.}
+    \label{fig:narrow_sweep}
+\end{figure}
+
+\begin{figure}[H]
+    \centering
+    \includegraphics[width=\linewidth]{images/particles_left_narrow_sweep_interactions.pdf}
+    \caption{Exploring the behavior of particles, where the amplitude of time-dependent potential $f = [0.1, 0.4, 0.7]$, for angular frequency $\omega_{V} \in (1.1, 1.7)$ MHz with particle interactions.}
+    \label{fig:narrow_sweep_interactions}
 \end{figure}
 \end{document}