\documentclass[../main.tex]{subfiles}
\resources % search for graphics in ../res/imgs/


\begin{document}

\section{Methods}

Newton's second law of motion is defined as
\begin{equation}
    \frac{d^2\vb{r}}{dt^2} = \frac{\vb{F}\left(t, \frac{d\vb{r}}{dt}, \vb{r}\right)}{m}.
    \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 redefine Newton's second law of motion as
\begin{equation}
    \begin{split}
        \frac{d\vb{r}}{dt} &= \vb{v} \\
        \frac{d\vb{v}}{dt} &= \frac{\vb{F}(t, \vb{v}, \vb{r})}{m}.
    \end{split}
    \label{eq:coupled}
\end{equation}

\subsection{Forward Euler}

For a single particle, the forward Euler method for a coupled system is 
expressed as

\begin{equation}
    \begin{split}
        \vb{r}_{i+1} &= \vb{r}_i + h \cdot \frac{d\vb{r}_i}{dt} = \vb{r}_i + h \cdot \vb{v}_i \\
        \vb{v}_{i+1} &= \vb{v}_i + h \cdot \frac{\vb{v}_i}{dt} = \vb{v}_i + h \cdot \frac{\vb{F}\left( t_i, \vb{v}_i, \vb{r}_i \right)}{m},
    \end{split}
\end{equation}
where $\vb{v}_i$ is the current velocity of the particle, $\vb{r}_i$ is the
current position of the particle, $m$ is the mass of the particle, and $h$
is a predetermined timestep.

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}
        \begin{algorithmic}
            \Procedure{Evolve forward Euler}{$particles, dt$}
            \State $new\_state \leftarrow particles$ \Comment{Create a copy of the particles}
            \For{ $i = 1, 2, \ldots , \left| particles \right|$ }
            \State $new\_state_i.\vb{v} \leftarrow particles_i.\vb{v} + dt \cdot \frac{\vb{F}}{m}$
            \State $new\_state_i.\vb{r} \leftarrow particles_i.\vb{r} + dt \cdot particles_i.\vb{v}$
            \EndFor
            \State $particles \leftarrow new\_state$
            \EndProcedure
        \end{algorithmic}
    \end{algorithm}
\end{figure}

\subsection{4th order Runge-Kutta}

For a single particle, we can express the 4th order Runge-Kutta (RK4) method as
\begin{equation}
    \begin{split}
        \vb{v}_{i+1} &= \vb{v}_i + \frac{h}{6} \left( \vb{k}_{\vb{v},1} 
            + 2\vb{k}_{\vb{v},2} + 2\vb{k}_{\vb{v},3} + \vb{k}_{\vb{v},4} 
            \right) \\
        \vb{r}_{i+1} &= \vb{r}_i + \frac{h}{6} \left( \vb{k}_{\vb{r},1} 
            + 2\vb{k}_{\vb{r},2} + 2\vb{k}_{\vb{r},3} + \vb{k}_{\vb{r},4} 
            \right),
    \end{split}
\end{equation}
where
\begin{equation}
    \begin{split}
    \vb{k}_{\vb{v},1} &= \frac{\vb{F}(t_i, \vb{v}_i, \vb{r}_i)}{m} \\
    \vb{k}_{\vb{r},1} &= \vb{v}_i \\
    \vb{k}_{\vb{v},2} &= \frac{\vb{F}(t_i+\frac{h}{2}, \vb{v}_i + h \cdot \frac{\vb{k}_{\vb{v},1}}{2}, \vb{r}_i + h \cdot \frac{\vb{k}_{\vb{r},1}}{2})}{m} \\
    \vb{k}_{\vb{r},2} &= \vb{v}_i + h \cdot \frac{\vb{k}_{\vb{v},1}}{2} \\
    \vb{k}_{\vb{v},3} &= \frac{\vb{F}(t_i+\frac{h}{2}, \vb{v}_i + h \cdot \frac{\vb{k}_{\vb{v},2}}{2}, \vb{r}_i + h \frac{\cdot \vb{k}_{\vb{r},2}}{2})}{m} \\
    \vb{k}_{\vb{r},3} &= \vb{v}_i + h \cdot \frac{\vb{k}_{\vb{v},2}}{2} \\
    \vb{k}_{\vb{v},4} &= \frac{\vb{F}(t_i+h, \vb{v}_i + h \cdot \vb{k}_{\vb{v},3}, \vb{r}_i + h \cdot \vb{k}_{\vb{r},3})}{m} \\
    \vb{k}_{\vb{r},4} &= \vb{v}_i + h \cdot \frac{\vb{k}_{\vb{v},1}}{2}.
    \end{split}
\end{equation}

\newpage

\begin{figure}
    \begin{algorithm}[H]
    \caption{RK4 method}
    \label{algo:rk4}
        \begin{algorithmic}
            \Procedure{Evolve RK4}{$dt$}
            \State $N \leftarrow \text{Number of particles inside the Penning trap}$
            \State $orig\_p \leftarrow \text{Copy of particles}$ \Comment{Keep the original particles}
            \State $tmp\_p \leftarrow \text{Array of particles of size }N$
            \State $\vb{k}_{\vb{r}} \leftarrow \text{2D array of vectors of size } 4 \cross N$
            \State $\vb{k}_{\vb{v}} \leftarrow \text{2D array of vectors of size } 4 \cross N$

            \For{ $i = 1, 2, \ldots, N$ }
            \State $\vb{k}_{\vb{r},1,i} \leftarrow \vb{v}$
            \State $\vb{k}_{\vb{v},1,i} \leftarrow \frac{\vb{F}}{m}$

            \State $tmp\_p_i.\vb{r} \leftarrow orig\_p_i.\vb{r} + \frac{dt}{2} \cdot \vb{k}_{\vb{r},1,i}$
            \State $tmp\_p_i.\vb{v} \leftarrow orig\_p_i.\vb{v} + \frac{dt}{2} \cdot \vb{k}_{\vb{v},1,i}$
            \EndFor

            \State $particles \leftarrow tmp\_p$ \Comment{Update particles}

            \For{ $i = 1, 2, \ldots, N$ }
            \State $\vb{k}_{\vb{r},2,i} \leftarrow \vb{v}$
            \State $\vb{k}_{\vb{v},2,i} \leftarrow \frac{\vb{F}}{m}$

            \State $tmp\_p_i.\vb{r} \leftarrow orig\_p_i.\vb{r} + \frac{dt}{2} \cdot \vb{k}_{\vb{r},2,i}$
            \State $tmp\_p_i.\vb{v} \leftarrow orig\_p_i.\vb{v} + \frac{dt}{2} \cdot \vb{k}_{\vb{v},2,i}$
            \EndFor

            \State $particles \leftarrow tmp\_p$ \Comment{Update particles}

            \For{ $i = 1, 2, \ldots, N$ }
            \State $\vb{k}_{\vb{r},3,i} \leftarrow \vb{v}$
            \State $\vb{k}_{\vb{v},3,i} \leftarrow \frac{\vb{F}}{m}$

            \State $tmp\_p_i.\vb{r} \leftarrow orig\_p_i.\vb{r} + dt \cdot \vb{k}_{\vb{r},3,i}$
            \State $tmp\_p_i.\vb{v} \leftarrow orig\_p_i.\vb{v} + dt \cdot \vb{k}_{\vb{v},3,i}$
            \EndFor

            \State $particles \leftarrow tmp\_p$ \Comment{Update particles}

            \For{ $i = 1, 2, \ldots, N$ }
            \State $\vb{k}_{\vb{r},4,i} \leftarrow \vb{v}$
            \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} 
            + \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} 
            + \vb{k}_{\vb{v},3,i} + \vb{k}_{\vb{v},4,i} \right)$
            
            \EndFor

            \State $particles \leftarrow tmp\_p$ \Comment{Final update}
            \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}


\biblio
\end{document}
