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


\begin{document}

\section{Methods}

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}

Newton's second law for a particle $i$ is then
\begin{equation}
    \frac{d^2\vb{r}_i}{dt^2} = \frac{\vb{F}\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}(t, \vb{v}_i, \vb{r}_i)}{m_i}.
    \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 $i$ is the current time step of the particle $m$ is the mass of the particle, and $h$
is the time 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}
        \begin{algorithmic}
            \Procedure{Evolve forward Euler}{$particles, dt$}
            \State $N \leftarrow \text{Number of particles in } particles$
            \State $a \leftarrow \text{Calculate } \frac{\vb{F_i}}{m_i} \text{ for each particle in } particles$
            \For{ $i = 1, 2, \ldots , N$ }
            \State $particles_i.\vb{r} \leftarrow particles_i.\vb{r} + dt \cdot particles_i.\vb{v}$
            \State $particles_i.\vb{v} \leftarrow particles_i.\vb{v} + dt \cdot a_i$
            \EndFor
            \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}

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.
around 8 separate loops to calculate all $\vb{k}$ values

\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}