Started the rewrite of the manual to correspond with the refactor of the code. Stream 1 is completely written into the new structure, now 9 other streams are still left to do

tex-docs/CLAuDE.tex

@@ -35,16 +35,12 @@ will be presented in SI units \cite{SI} between brackets like: $T$: The temperat
 to relate SI units to your preferred system of units, please refer to the internet for help with that. There are great calculators online where you only need to plug in a number and select the 
 right units.
 
-Within this manual we will not concern ourselves with plotting the data, instead we focus on the physics side of things and translating formular into code. If you are interested in how the
+Within this manual we will not concern ourselves with plotting the data, instead we focus on the physics side of things and translating formulae into code. If you are interested in how the
 plotting of the data works, or how loading and saving data works, please refer to the relevant stream on Simon's Twitch page \cite{twitch}.
 
-This manual is for the toy model, which is as of now still in development. Therefore this manual will be in chronological order, explaining everything in the same order as it has been done. 
+This manual is for the toy model, which is as of now still in development. One important thing to note is that the layout may change significantly when new sections are added. This is due to the amount of code that is added/changed. If a lot of code changes, a lot of so called
-There are plans to eventually modularise the whole model into seperate parts that can be extended by the community. When that will hit development, a new manual for that version of the model
-will be made that treats things per topic instead of chronological order. 
-
-One important thing to note, the layout may change significantly when new sections are added. This is due to the amount of code that is added/changed. If a lot of code changes, a lot of so called
 'algorithm' blocks are present which have different placement rules than just plain text. Therefore it may occur that an algorithm is referenced even though it is one or two pages later. This is 
-a pain to fix and if something later on changes the whole layout may be messed up again and is a pain to fix again. Hence I opt to let \LaTeX (the software/typeset language used to create this 
+a pain to fix and if something later on changes, the whole layout may be messed up again and is a pain to fix again. Hence I opt to let \LaTeX (the software/typeset language used to create this 
 manual) figure out the placement of the algorithm blocks, which may or may not be in the right places.
 
 Lastly, the manual is now up on the Planet Factory GitHub repository\cite{claudeGit}, together with all the source code. There is also a fork \cite{nomGit} that also contains the source code. 
@@ -53,8 +49,6 @@ particular stream is missing in the version on the Planet Factory repository, ch
 patient, or you can start writing a part of the manual yourself! Don't forget to ping me in the Discord to notify me of any additions (GitHub refuses to send me emails so I have no other way of 
 knowing).
 
-\input{streams/Stream1.tex}
-
 \input{streams/Stream2.tex}
 
 \input{streams/Stream3.tex}
@@ -73,6 +67,16 @@ knowing).
 
 \input{streams/Stream10.tex}
 
+\input{topics/control_panel.tex}
+
+\input{topics/util_funcs.tex}
+
+\input{topics/radiation.tex}
+
+%Velocity
+
+%Advection
+
 \newpage
 \input{streams/TTNMETAF.tex}
 

tex-docs/references.bib

@@ -45,6 +45,17 @@ pages={1328},
 edition={14th global}
 }
 
+@inbook{specificHeat, 
+place={Harlow}, 
+title={Sears and Zemanskys University physics with modern physics}, 
+publisher={Pearson Education}, 
+author={Young, Hugh D. and Freedman, Roger A. and Ford, A. Lewis}, 
+year={2016},
+chapter={17},
+pages={581},
+edition={14th global}
+}
+
 @inbook{thermo1, 
 place={Harlow}, 
 title={Sears and Zemanskys University physics with modern physics}, 

tex-docs/streams/Stream2.tex

@@ -3,7 +3,7 @@ In its current state, CLaUDE has a static planet. This means that the planet rem
 themselves. But before we start adding layers, let's talk about a term you will hear more often: numerical instability.
 
 Numerical instability occurs when you first run the model. This is due to the nature of the equations. Nearly all equations are continuous, which means that they are always at work. However 
-when you start the model, the equations were not at work yet. It is as if you suddenly give a random meteor an atmosphere, place it in orbit around a star and don't touch i for a bit. You will 
+when you start the model, the equations were not at work yet. It is as if you suddenly give a random meteor an atmosphere, place it in orbit around a star and don't touch it for a bit. You will 
 see that the whole system oscilates wildly as it adjusts to the sudden changes and eventually it will stabilise. Another term you might encounter is blow up, this occurs when when the model 
 suddenly no longer behaves like it should. This is most likely caused by mistakes in the code or incorrect paramter initialisation. Be wary of the existence of both factors, and do not dismiss 
 a model if it behaves weirdly as it has just started up.

tex-docs/streams/Stream3.tex

@@ -46,35 +46,6 @@ $u$ in both $x$ and $y$ directions. Then if we write out $\nabla u$ we get \auto
     \end{equation}
 \end{subequations}
 
-Now that we have the momentum equations sorted out, we need to define a method to do the gradient calculations for us. Therefore we define two functions \autoref{alg:gradient x} and 
-\autoref{alg:gradient y} that calculate the $x$ and $y$ gradients respectively.
-
-\begin{algorithm}[hbt]
-    \SetKwInOut{Input}{Input}
-    \SetKwInOut{Output}{Output}
-    \Input{Matrix (double array) $a$, first index $i$, second index $j$}
-    \Output{Gradient in the $x$ direction}
-    $grad \leftarrow \frac{a[i, (j + 1)\text{ mod } nlon] - a[i, (j - 1) \text{ mod } nlon]}{\delta x[i]}$ \;
-    \Return{$grad$} \;
-    \caption{Calculating the gradient in the $x$ direction}
-    \label{alg:gradient x}
-\end{algorithm}
-
-\begin{algorithm}[hbt]
-    \SetKwInOut{Input}{Input}
-    \SetKwInOut{Output}{Output}
-    \Input{Matrix (double array) $a$, first index $i$, second index $j$}
-    \Output{Gradient in the $y$ direction}
-    \eIf{$i == 0$ or $i == nlat - 1$}{
-        $grad \leftarrow 0$ \;
-    }{
-        $grad \leftarrow \frac{a[i + 1, j] - a[i - 1 j]}{\delta y}$ \;
-    }
-    \Return $grad$ \;
-    \caption{Calculating the gradient in the $y$ direction}
-    \label{alg:gradient y}
-\end{algorithm}
-
 With the gradient functions defined, we can move on to the main code for the momentum equations. The main loop is shown in \autoref{alg:stream3}. Do note that this loop replaces the one
 in \autoref{alg:stream2v2} as these calculate the same thing, but the new algorithm does it better.
 
@@ -117,8 +88,7 @@ The diffusion equation, as written in \autoref{eq:diffusion}, describes how the
 \end{equation}
 
 Now to get this into code we need the following algorithms \autoref{alg:laplacian} and \autoref{alg:diffusion}. \autoref{alg:laplacian} implements the laplacian operator, whereas 
-\autoref{alg:diffusion} implements the diffusion calculations. $\Delta_x$ and $\Delta_y$ in \autoref{alg:laplacian} represents the calls to \autoref{alg:gradient x} and \autoref{alg:gradient y} 
+\autoref{alg:diffusion} implements the diffusion calculations.  $\nabla^2$ in \autoref{alg:diffusion} represents the call to \autoref{alg:laplacian}.
-respectively. $\nabla^2$ in \autoref{alg:diffusion} represents the call to \autoref{alg:laplacian}.
 
 \begin{algorithm}
     \SetKwInOut{Input}{Input}
@@ -164,10 +134,6 @@ The advection equation is shown in \autoref{eq:advection}. The symbols are:
     \label{eq:advection}
 \end{equation}
 
-As we expect to use the divergence operator more often throughout our model, let us define a seperate function for it in \autoref{alg:divergence}. $\Delta_x$ and $\Delta_y$ in 
-\autoref{alg:divergence} represents the calls to \autoref{alg:gradient x} and \autoref{alg:gradient y} respectively. We do the multiplication with the velocity vector here already, as we expect 
-that we might use it in combination with the divergence operator more frequently.
-
 \begin{algorithm}
     \SetKwInOut{Input}{Input}
     \SetKwInOut{Output}{Output}

tex-docs/streams/Stream5.tex

@@ -16,90 +16,7 @@ are no longer indexed by $lat, lon$ but are indexed by $lat, lon, layer$.
     \label{alg:more layers}
 \end{algorithm}
 
-We also need to change all the gradient functions (\autoref{alg:gradient x} and \autoref{alg:gradient y}) to incorporate the atmospheric layers. Additionally we need a new gradient function that 
-calculates the gradient in the $z$ direction (vertical). Let us first change the existing gradient functions to take the atmospheric layer in effect. The changes can be found in 
-\autoref{alg:gradient x layer} and \autoref{alg:gradient y layer}. Let us improve the gradient in the $y$ direction as well. Since we are using the central difference method (calculating the
-gradient by taking the difference of the next grid cell and the previous grid cell) there is no gradient at the poles. What we can do instead of returning $0$ for those cases is forward 
-differencing (calculating the gradient by taking the difference of the cell and the next/previous cell, multiplied by $2$ to keep it fair). A special change in both functions is checking whether
-$k$ is equal to \texttt{NULL}. We do this as sometimes we want to use this function for matrices that does not have the layer dimension. Hence we define a default value for $k$ which is 
-\texttt{NULL}. \texttt{NULL} is a special value in computer science. It represents nothing. This can be useful sometimes if you declare a variable to be something but it is referring to something
-that has been deleted or it is returned when some function fails. It usually indicates that something special is going on. So here we use it in the special case where we do not want to consider
-the layer part in the gradient.
 
-\begin{algorithm}[hbt]
-    \SetKwInOut{Input}{Input}
-    \SetKwInOut{Output}{Output}
-    \Input{Matrix (double array) $a$, first index $i$, second index $j$, third index $k$ with default value \texttt{NULL}}
-    \Output{Gradient in the $x$ direction}
-    \eIf{$k == \texttt{NULL}$}{
-        $grad \leftarrow \frac{a[i, (j + 1)\text{ mod } nlon] - a[i, (j - 1) \text{ mod } nlon]}{\delta x[i]}$ \;
-    }{
-        $grad \leftarrow \frac{a[i, (j + 1)\text{ mod } nlon, k] - a[i, (j - 1) \text{ mod } nlon, k]}{\delta x[i]}$ \;
-    }
-    \Return{$grad$} \;
-    \caption{Calculating the gradient in the $x$ direction}
-    \label{alg:gradient x layer}
-\end{algorithm}
-
-\begin{algorithm}[hbt]
-    \SetKwInOut{Input}{Input}
-    \SetKwInOut{Output}{Output}
-    \Input{Matrix (double array) $a$, first index $i$, second index $j$, third index $k$ with default value \texttt{NULL}}
-    \Output{Gradient in the $y$ direction}
-    \eIf{$k == \texttt{NULL}$}{
-        \uIf{$i == 0$}{
-            $grad \leftarrow 2 \frac{a[i + 1, j] - a[i, j]}{\delta y}$ \;
-        }\uElseIf{$i == nlat - 1$}{
-            $grad \leftarrow 2 \frac{a[i, j] - a[i - 1, j]}{\delta y}$ \;
-        }\uElse{
-            $grad \leftarrow \frac{a[i + 1, j] - a[i - 1 j]}{\delta y}$ \;
-        }
-    }{
-        \uIf{$i == 0$}{
-            $grad \leftarrow 2 \frac{a[i + 1, j, k] - a[i, j, k]}{\delta y}$ \;
-        }\uElseIf{$i == nlat - 1$}{
-            $grad \leftarrow 2 \frac{a[i, j, k] - a[i - 1, j, k]}{\delta y}$ \;
-        }\uElse{
-            $grad \leftarrow \frac{a[i + 1, j] - a[i - 1 j]}{\delta y}$ \;
-        }
-    }
-    \Return $grad$ \;
-    \caption{Calculating the gradient in the $y$ direction}
-    \label{alg:gradient y layer}
-\end{algorithm}
-
-With those changes done, let us define the gradient in the $z$ direction. The function can be found in \autoref{alg:gradient z layer}. Here $a.dimensions$ is the attribute that tells us how 
-deeply nested the array $a$ is. If the result is $1$ we have just a normal array, if it is $2$ we have a double array (an array at each index of the array) which is also called a matrix and if it 
-is $3$ we have a triple array. We need this because we have a one-dimensional case, for when we do not use multiple layers and a three-dimensional case for when we do use multiple layers. This 
-distinction is needed to avoid errors being thrown when running the model with one or multiple layers.
-
-\begin{algorithm}[hbt]
-    \SetKwInOut{Input}{Input}
-    \SetKwInOut{Output}{Output}
-    \Input{Matrix (double array) $a$, first index $i$, second index $j$, third index $k$}
-    \Output{Gradient in the $z$ direction}
-    \uIf{$a.dimensions == 1$}{
-        \uIf{$k == 0$}{
-            $grad \leftarrow \frac{a[k + 1] - a[k]}{\delta z[k]}$ \;
-        }\uElseIf{$k == nlevels - 1$}{
-            $grad \leftarrow \frac{a[k] - a[k - 1]}{\delta z[k]}$ \;
-        }\uElse{
-            $grad \leftarrow \frac{a[k + 1] - a[k - 1]}{2\delta z[k]}$ \;
-        }
-    } \uElse {
-        \uIf{$k == 0$}{
-            $grad \leftarrow \frac{a[i, j, k + 1] - a[i, j, k]}{\delta z[k]}$ \;
-        }\uElseIf{$k == nlevels - 1$}{
-            $grad \leftarrow \frac{a[i, j, k] - a[i, j, k - 1]}{\delta z[k]}$ \;
-        }\uElse{
-            $grad \leftarrow \frac{a[i, j, k + 1] - a[i, j, k - 1]}{2\delta z[k]}$ \;
-        }
-    }
-    
-    \Return $grad$ \;
-    \caption{Calculating the gradient in the $z$ direction}
-    \label{alg:gradient z layer}
-\end{algorithm}
 
 As you can see, we have used $\delta z$ however, we have not defined it yet. Let us do that in \autoref{alg:gradient z}.
 
@@ -112,58 +29,10 @@ As you can see, we have used $\delta z$ however, we have not defined it yet. Let
     \label{alg:gradient z}
 \end{algorithm}
 
-Let's incorporate the changes for the Laplacian operator (\autoref{alg:laplacian}) as well. The new code can be found in \autoref{alg:laplacian layer}.
-
-\begin{algorithm}[hbt]
-    \SetKwInOut{Input}{Input}
-    \SetKwInOut{Output}{Output}
-    \Input{A matrix (double array) a}
-    \Output{A matrix (double array) with results for the laplacian operator for each element}
-    \eIf{$a.dimensions == 2$}{
-        \For{$lat \in [1, nlat - 1]$}{
-            \For{$lon \in [0, nlon]$}{
-                $output[lat, lon] \leftarrow \frac{\Delta_x(a, lat, (lon + 1) \text{ mod } nlon) - \Delta_x(a, lat, (lon - 1) \text{ mod } nlon)}{\delta x[lat]} + \frac{\Delta_y(a, lat + 1, lon) - 
-                \Delta_y(a, lat - 1, lon)}{\delta y}$\;
-            }
-        }
-    }{
-        \For{$lat \in [1, nlat - 1]$}{
-            \For{$lon \in [0, nlon]$}{
-                \For{$k \in [0, nlevels - 1]$}{
-                    $output[lat, lon, k] \leftarrow \frac{\Delta_x(a, lat, (lon + 1) \text{ mod } nlon, k) - \Delta_x(a, lat, (lon - 1) \text{ mod } nlon, k)}{\delta x[lat]} + \frac{\Delta_y(a, 
-                    lat + 1, lon, k) - \Delta_y(a, lat - 1, lon, k)}{\delta y} + \frac{\Delta_z(a, lat, lon, k + 1) - \Delta_z(a, lat, lon, k + 1)}{2\delta z[k]}$\;
-                }
-            }
-        }
-    }
-    
-    \Return{$ouput$} \;
-    \caption{Calculate the laplacian operator over a matrix a}
-    \label{alg:laplacian layer}
-\end{algorithm}
-
 Of course we also need to incorporate the new layers in the divergence operator (\autoref{alg:divergence}). The new changes can be found in \autoref{alg:divergence layer}. Here we use $w$, the 
 vertical wind velocity. We define $w$ in the same way as $u$ and $v$, it is all zeroes (in the beginning) and has the same dimensions as $u$ and $v$.
 
-\begin{algorithm}[!hbt]
+
-    \SetKwInOut{Input}{Input}
-    \SetKwInOut{Output}{Output}
-    \Input{A matrix (double array) $a$}
-    \Output{A matrix (double array) containing the result of the divergence operator taken over that element}
-    $dim_1 \leftarrow \text{ Length of } a \text{ in the first dimension}$ \;
-    \For{$i \in [0, dim_1]$}{
-        $dim_2 \leftarrow \text{ Length of } a \text{ in the second dimension (i.e. the length of the array stored at index } i)$ \;
-        \For{$j \in [0, dim_2]$}{
-            $dim_3 \leftarrow \text{ Length of } a \text{ in the third dimension}$ \;
-            \For{$k \in [0, dim_3]$}{
-                $output[i, j] \leftarrow \Delta_x(au, i, j, k) + \Delta_y(av, i, j, k) + \Delta_z(aw, i, j, k)$ \;
-            }
-        }
-    }
-    \Return{$output$} \;
-    \caption{Calculate the result of the divergence operator on a vector}
-    \label{alg:divergence layer}
-\end{algorithm}
 
 With all those changes in the functions done, let us incorporate the changes into the model itself. We now need to account for the temperature change throughout the layers. Let us look at the 
 atmospheric temperature equation again (\autoref{eq:atmos change}). We need to account for one more thing, the absorbtion of energy from another layer. The new equation is shown in 

tex-docs/streams/TTNMETAF.tex

@@ -5,51 +5,4 @@
 Potential is the energy change that occurs when the position of an object changes \cite{potential}. There are many potentials, like electric potential, gravitational potential and elastic 
 potential. Let me explain the concept with an example. Say you are walking on a set of stairs in the upwards direction. As your muscles move to bring you one step upwards, energy that is used
 by your muscles is converted into gravitational potential. Now imagine you turn around and go downwards instead. Notice how that is easier? That is due to the gravitational potential being 
-converted back into energy so your muscles have to deliver less energy to get you down. The potential is usually tied to a force, like the gravitational force.
+converted back into energy so your muscles have to deliver less energy to get you down. The potential is usually tied to a force, like the gravitational force.
-
-\subsection{Laplacian Operator} \label{sec:laplace}
-The Laplacian operator ($\nabla^2$, sometimes also seen as $\Delta$) has two definitions, one for a vector field and one for a scalar field. The two concepts are not indpendent, a vector field 
-is composed of scalar fields \cite{vectorscalarfields}. Let us define a vector field first. A vector field is a function whose domain and range are a subset of the Eucledian $\mathbb{R}^3$ space. 
-A scalar field is then a function consisting out of several real variables (meaning that the variables can only take real numbers as valid values). So for instance the circle equation 
-$x^2 + y^2 = r^2$ is a scalar field as $x, y$ and $r$ are only allowed to take real numbers as their values. 
-
-With the vector and scalar fields defined, let us take a look at the Laplacian operator. For a scalar field $\phi$ the laplacian operator is defined as the divergence of the gradient of $\phi$
-\cite{laplacian}. But what are the divergence and gradient? The gradient is defined in \autoref{eq:gradient} and the divergence is defined in \autoref{eq:divergence}. Here $\phi$ is a vector 
-with components $x, y, z$ and $\Phi$ is a vector field with components $x, y, z$. $\Phi_1, \Phi_2$ and $\Phi_3$ refer to the functions that result in the corresponding $x, y$ and $z$ values 
-\cite{vectorscalarfields}. Also, $i, j$ and $k$ are the basis vectors of $\mathbb{R^3}$, and the multiplication of each term with their basis vector results in $\Phi_1, \Phi_2$ and $\Phi_3$
-respectively. If we then combine the two we get the Laplacian operator, as in \autoref{eq:laplacian scalar}.
-
-\begin{subequations}
-    \begin{equation}
-        \text{grad } \phi = \nabla \phi = \frac{\delta \phi}{\delta x}i + \frac{\delta \phi}{\delta y}j + \frac{\delta \phi}{\delta z}k
-        \label{eq:gradient}
-    \end{equation}
-    \begin{equation}
-        \text{div} \Phi = \nabla \cdot \Phi = \frac{\delta \Phi_1}{\delta x} + \frac{\delta \Phi_2}{\delta y} + \frac{\delta \Phi_3}{\delta z}
-        \label{eq:divergence}
-    \end{equation}
-    \begin{equation}
-        \nabla^2 \phi = \nabla \cdot \nabla \phi = \frac{\delta^2 \phi}{\delta x^2} + \frac{\delta^2 \phi}{\delta y^2} + \frac{\delta^2 \phi}{\delta z^2}
-        \label{eq:laplacian scalar}
-    \end{equation}
-\end{subequations}
-
-For a vector field $\Phi$ the Laplacian operator is defined as in \autoref{eq:laplacian vector}. Which essential boils down to taking the Laplacian operator of each function and multiply it by
-the basis vector.
-
-\begin{equation}
-    \nabla^2 \Phi = (\nabla^2 \Phi_1)i + (\nabla^2 \Phi_2)j + (\nabla^2 \Phi_3)k
-    \label{eq:laplacian vector}
-\end{equation}
-
-\subsection{Interpolation} \label{sec:interpolation}
-Interpolation is a form of estimation, where one has a set of data points and desires to know the values of other data points that are not in the original set of data points\cite{interpolation}. 
-Based on the original data points, it is estimated what the values of the new data points will be. There are various forms of interpolation like linear interpolation, polynomial interpolation 
-and spline interpolation. The CLAuDE model uses linear interpolation which is specified in \autoref{eq:interpolation}. Here $z$ is the point inbetween the known data points $x$ and $y$. 
-$\lambda$ is the factor that tells us how close $z$ is to $y$ in the interval $[0, 1]$. If $z$ is very close to $y$, $\lambda$ will have the value on the larger end of the interval, like 0.9.
-Whereas if $z$ is close to $x$ then $\lambda$ will have a value on the lower end of the interval, like 0.1.
-
-\begin{equation}
-    z = (1 - \lambda)x + \lambda y
-    \label{eq:interpolation}
-\end{equation}

tex-docs/topics/control_panel.tex

@@ -0,0 +1,104 @@
+\section{Control Panel} \label{sec:cp}
+Before we dive in an start modelling the planet, let us first set up a control panel that will influence how the model will behave and effectively decides what type of planet we model.
+
+\subsection{The Beginning}
+In the beginning there was nothing, and then there was "Hello World!" Or at least that is how many projects start. Why? you might ask, which is a perfectly valid question. In Computer Science, 
+"Hello World!" is very simple code that we use to test whether all the tools we need to get coding works. This checks whether the computer compiles the code and is able to execute it and whether 
+the code editor (IDE, Integrated Development Environment) starts the right processes to get the code compiled and executed. Oh right we were talking about CLAuDE, ahem.
+
+Every project must have its beginning. And with CLAuDE I made the decision to start explaining the Control Panel first. This is to get you familiar with notation and to lay down some basics. To 
+do that we start with the fixed part of the Control Panel, the physical constants. Many things vary from planet to planet, how much radiation they receive from their star, how strong their 
+gravity is, how fast they spin around their axis and many many more. What does not change are the physical constants, well because they are constant. The Stefan-Boltzmann constant for instance
+does not change. Whether you are on Earth, in space or on Jupiter, the value of the Stefan-Boltzmann constant will remain the same. 
+
+The Stefan-Boltzmann constant is denoted by $\sigma$ and has a value of $5.670373 \cdot 10^-8 \ (Wm^{-2}K^{-4})$ \cite{stefan-boltzmann}. The $\sigma$ is a greek letter called sigma. Greek 
+letters are often used in mathematics, as well as in physics or any other discipline that relies on maths (spoiler alert, quite a lot). Treat it like a normal letter in maths, representing a 
+number that you either do not know yet or is too long or cumbersome to write down every time. The Stefan-Boltzmann constant is denoted in scientific notation, a number followed by the order of 
+magnitude. It is denoted as a multiplication, because that is what you have to do to get the real number. An example: $4.3 \cdot 10^2 = 430$ and $4.3 \cdot 10^{-2} = 0.043$. The
+letters behind the numbers are units, how we give meaning to the numbers. If I say that I am $1.67$ does not mean anything. Do I mean inches, centimeters, meters, miles? That is why we need units
+as they give meaning to the number. they tell us whether the number is talking about speed, distance, time, energy and many other things. In this manual we will use SI units. Behind all the
+letters you will find the following: [number]. This is a citation, a reference to an external source where you can check whether I can still read. If I pull a value out somewhere I will insert a 
+citation to show that I am not making these numbers up. This is what scientists use to back up their claims if they do not want to redo the work that others have done. I mean what is the point 
+of re-inventing the wheel if there is a tyre company next door? That is why scientists citate.
+
+So with that out of the way, let us write down some constants. Why do I do this here? Because a lot of constants are used everywhere and I am too lazy to relicate them every time. If you see a 
+letter or symbol that is not explicitly explained, then it is most likely a constant that we discuss here in the control panel. 
+
+\subsection{Physical Constants}
+As mentioned before, physical constants do not change based on where you are in the universe. Below you will find an overview of all the relevant constants together with their units. And a short
+explanation where they are used or what they represent. To see them in action, consult the other sections of this manual, you will find them in equations all throughout this document.
+
+\subsubsection{The Gas Constant}
+The Gas constant, $R = 8.3144621$ ($J(mol)^{-1}K$) \cite{idealGas} is the constant used to relate the temperature of the gas to the pressure and the volume. One would expect this constant to be 
+different per gas, but under high enough temperatures and low enough pressure the gas constant is the same for all gases. 
+
+\subsubsection{The Specific Heat Capacity}
+The specific heat capacity $c$ depicts how much energy is required to heat the object by one degree Kelvin per unit mass ($\frac{J}{Kg \cdot K}$) \cite{specificHeat}. This varies per material 
+and is usually indicated by a subscript. The specific heat capacity for water for instance is $c_w = 4190  JKg^{-1}K^{-1}$. Specific heat capacities also exist in the form of $Jg^{-1}K^{-1}$, 
+$Jmol^{-1}K^{-1}$ and $Jcm^{-3}K^{-1}$ which you can use in various circumstances, depending on what information you have.
+
+\subsubsection{Mole}
+Mole is the amount of particles ($6.02214076 \cdot 10^{23}$) in a substance, where the average weight of one mole of particles in grams is about the same as the weight of one particle in atomic 
+mass units ($u$)\cite{mole}. This is not a physical constant perse, but more like a unit ($mol$). Though it is still important enough to be added here for future reference. All other units are 
+way more intuitive and are assumed to be known.
+
+\subsubsection{The Stefan-Boltzmann Constant}
+The Stefan-Boltzmann constant, $\sigma = 5.670373 \cdot 10^-8 \ (Wm^{-2}K^{-4})$ \cite{stefan-boltzmann} is used in the Stefan-Boltzmann law (more on that in %insert reference here). 
+
+\subsection{Planet Specific Variables}
+The following set of variables vary per planet, that's why we call them variables since they vary. Makes sense right? We add them here as we will use them throughout the manual. The advantage
+of that is quite significant. If you want to test things for a different planet, you only need to change the values in one place, instead of all places where you use it. If there is one thing 
+that we computer scientists hate is doing work, we like being lazy and defining things in one place means that we can be lazy if we need to change it. So we put in the extra work now, so we do 
+not have to do the extra work in the future. That's actually a quite accurate description of computer scientists, doing hard work so that they can be lazy in the future.
+
+\subsubsection{The Passage of Time}
+On Earth we have various indications of how much time has passed. While most of them remain the same throughout the universe, like seconds, minutes and hours, others vary throughout the universe,
+like days, months and years. Here we specify how long the variable quantities of time are for the planet we want to consider as they are used in the code. This can be seen in 
+\autoref{alg:time constants}. Here a $\leftarrow$ indicates that we assign a value to the variable name before it, so that we can use the variable name in the code instead of the value, which 
+has the advantage I indicated before. // means that we start a comment, which is text that the code ignores and does not tell the cpu about. Not that the cpu would understand it, but that just 
+means less work for the computer. Yes computers are lazy too.
+
+\begin{algorithm*}
+    \caption{Definition of how much time it takes for a day and a year on a planet and how much time on the planet passes before we start another calculation run}
+    \label{alg:time constants}
+    \SetKwComment{Comment}{//}{}
+    $day \leftarrow 60*60*24$ \Comment*[l]{Length of one day in seconds ($s$)}
+    $year \leftarrow 365*day$ \Comment*[l]{Length of one year in seconds ($s$)}
+    $\delta t \leftarrow 60 * 9$ \Comment*[l]{How much time is between each calculation run in seconds ($s$)}
+\end{algorithm*}
+
+\subsubsection{The Planet Passport}
+Each planet is different, so why should they all have the same gravity? Oh wait, they don't. Just as they are not all the same size, tilted as much and their atmospheres differ. So here we define
+all the relevant variables that are unique to a planet, or well not necessarily unique but you get the idea. This can all be found in \autoref{alg:planet constants}.
+
+\begin{algorithm}
+    \caption{Defining the constants that are specific to a planet}
+    \label{alg:planet constants}
+    \SetKwComment{Comment}{//}{}
+    $g \leftarrow 9.81$ \Comment*[l]{Magnitude of gravity on the planet in $ms^{-2}$}
+    $\alpha \leftarrow -23.5$ \Comment*[l]{By how many degrees the planet is tilted with respect to the star's plane}
+    $top \leftarrow 50*10^3$ \Comment*[l]{How high the top of the atmosphere is with respect to the planet surface in meters ($m$)}
+    $ins \leftarrow 1370$ \Comment*[l]{Amount of energy from the star that reaches the planet per unit area ($Jm^{-2}$)}
+    $\epsilon \leftarrow 0.75$ \Comment*[l]{Absorbtivity of the atmosphere, fraction of how much of the total energy is absorbed (unitless)}
+    %$R \leftarrow 6.4*10^6$ \Comment*[l]{The radius of the planet in meters ($m$)}
+\end{algorithm}
+
+\subsubsection{Model Specific Parameters}
+These parameters cannot be found out in the wild, they only exist within our model. They control things like the size of a cell on the latitude longitude grid (more on that in later sections),
+how much time the model gets to spin up. We need the model to spin up in order to avoid numerical instability. Numerical instability occurs when you first run the model. This is due to the nature of the equations. Nearly all equations are continuous, which means that they are always at work. However 
+when you start the model, the equations were not at work yet. It is as if you suddenly give a random meteor an atmosphere, place it in orbit around a star and don't touch it for a bit. You will 
+see that the whole system oscilates wildly as it adjusts to the sudden changes and eventually it will stabilise. We define the amount of time it needs to stabilise as the spin up time. All 
+definitions can be found in \autoref{alg:model constants}. What the $adv$ boolean does is enabling or disabling advection, a process described in <section> which does not work yet.
+
+\begin{algorithm}
+    \caption{Defining the paramters that only apply to the model}
+    \label{alg:model constants}
+    \SetKwComment{Comment}{//}{}
+    $resolution \leftarrow 3$ \Comment*[l]{The amount of degrees on the latitude longitude grid that each cell has, with this setting each cell is 3 degrees latitude high and 3 degrees 
+    longitude wide}
+    $nlevels \leftarrow 10$ \Comment*[l]{The amount of layers in the atmosphere}
+    $\delta t_s \leftarrow 60*137$ \Comment*[l]{The time between calculation rounds during the spin up period in seconds ($s$)}
+    $t_s \leftarrow 5*day$ \Comment*[l]{How long we let the planet spin up in seconds ($s$)}
+    $adv \leftarrow \texttt{FALSE}$ \Comment*[l]{Whether we want to enable advection or not}
+    $adv\_boun \leftarrow 8$ \Comment*[l]{How many cells away from the poles where we want to stop calculating the effects of advection}
+\end{algorithm}

tex-docs/topics/radiation.tex

@@ -1,19 +1,18 @@
-\section{The Beginning}
+\section{Radiation}
 \subsection{The First Law of Thermodynamics and the Stefan-Boltzmann Equation}
-The beginning of CLAuDE is based upon one of the most important laws of physics: "Energy is neither created nor destroyed, only changed from one form to another." In otherwords, if energy goes into an object it must 
+The beginning of CLAuDE is based upon one of the most important laws of physics: "Energy is neither created nor destroyed, only changed from one form to another." In otherwords, if energy goes 
-equal the outflowing energy plus the change of internal energy. This is captured in Stefan-Boltzmann's law (\autoref{eq:stefan-boltzmann}) \cite{stefan-boltzmann}. 
+into an object it must equal the outflowing energy plus the change of internal energy. This is captured in Stefan-Boltzmann's law (\autoref{eq:stefan-boltzmann}) \cite{stefan-boltzmann}. 
 
 Here we assume that the planet is a black body, i.e. it absorbs all radiation (energy waves, some waves are visible like light, others are invisible like radio signals) on all wavelengths. 
 In \autoref{eq:stefan-boltzmann} the symbols are:
 
 \begin{itemize}
-    \item $S$: The energy that reaches the top of the atmosphere, coming from the sun or a similar star, per meters squared $Jm^{-2}$. This is also called the insolation.
+    \item $S$: The energy that reaches the top of the atmosphere, coming from the sun or a similar star, per second per meters squared $Wm^{-2}$. This is also called the insolation.
     \item $\sigma$: The Stefan-Boltzmann constant, $5.670373 \cdot 10^-8 \ (Wm^{-2}K^{-4})$ \cite{stefan-boltzmann}.
     \item $T$: The temperature of the planet ($K$).
 \end{itemize}
 
-Technically speaking \autoref{eq:stefan-boltzmann} is incorrect, as there is a mismatch in units. However, that is corrected in \autoref{eq:basis sphere final} so there is no need to worry about 
+The energy difference between the energy that reaches the atmosphere and the temperature of the planet must be equal to the change in temperature of the planet, which is written in 
-it just yet. The energy difference between the energu that reaches the atmosphere and the temperature of the planet must be equal to the change in temperature of the planet, which is written in 
 \autoref{eq:sb rewritten}. The symbols on the right hand side are:
 
 \begin{itemize}
@@ -22,8 +21,9 @@ it just yet. The energy difference between the energu that reaches the atmospher
     \item $\Delta T$: The change in temperature ($K$).
 \end{itemize}
 
-We want to know the change of temperature $\Delta T$, so we rewrite the equation into \autoref{eq:sb rewritten2}. Here we added the $\delta t$ term to account for the time difference (or time step). This is needed as
+We want to know the change of temperature $\Delta T$, so we rewrite the equation into \autoref{eq:sb rewritten2}. Here we added the $\delta t$ term to account for the time difference (or time 
-we need an interval to calculate the difference in temperature over. Also we needed to make the units match, and by adding this time step the units all match up perfectly.
+step). This is needed as we need an interval to calculate the difference in temperature over. Also we need to get the energy that we get ($J$) and not the energy per second ($W$), and by adding 
+this time step the units all match up perfectly.
 
 \begin{subequations}
     \begin{equation}
@@ -122,52 +122,50 @@ Pseudocode is a representation of real code. It is meant to be an abstraction of
 it in their language of preference. This is usually easier to read than normal code, but more difficult to read than mathematical formulae. If you are unfamiliar with code or coding, look up 
 a tutorial online as there are numerous great ones.
 
-The pseudocode in \autoref{alg:stream1v1} defines the main loop of the model. All values are initialised beforehand, based on either estimations, trial and error or because they are what they 
+The pseudocode in \autoref{alg:stream1v1} defines the main function of the radiation part of the model. All values are initialised beforehand, based on either estimations, trial and error or 
-are (like the Stefan-Boltzmann constant $\sigma$). The total time $t$ starts at 0 and increases by $\delta t$ after every update of the temperature. This is to account for the total time that 
+because they are what they are (like the Stefan-Boltzmann constant $\sigma$), which is all done in \autoref{sec:cp}. The total time $t$ starts at 0 and increases by $\delta t$ after every 
-the model has simulated (and it is also used later). What you may notice is the $T_p[lan, lon]$ notation. This is to indicate that $T_p$ saves a value for each $lan$ and $lon$ combination.
+update of the temperature. This is to account for the total time that the model has simulated (and it is also used later). What you may notice is the $T_p[lan, lon]$ notation. This is to indicate 
-It is initialised as all zeroes for each index pair, and the values is changed based on the calculations. You can view $T_p$ like the whole latitude longitude grid, where $T_p[lat, lon]$ is an 
+that $T_p$ saves a value for each $lan$ and $lon$ combination. It is initialised as all zeroes for each index pair, and the values is changed based on the calculations. You can view $T_p$ like 
-individual cell of that grid indexed by a specific latitude longitude combination.
+the whole latitude longitude grid, where $T_p[lat, lon]$ is an individual cell of that grid indexed by a specific latitude longitude combination. Do note that from here on most, if not all
+functions need to be called \footnote{In case you are unfamiliar with calls, defining a function is defining how it works and calling a function is actually using it.} from another file which I 
+will call the master. The master file decides what parts of the model to use, what information it uses for plots and the like. We do it this way because we want to be able to switch out 
+calculations. Say that I find a more efficient way, or more detailed way, to calculate the temperature change. If everything was in one file, then I need to edit the source code of the project.
+With the master file structure, I can just swap out the reference to the project's implementation with a reference to my own implementation. This makes the life of the user (in this case the 
+programmer who has another implementation) easier and makes changing calculations in the future easier as well. Also note that what we pass on as parameters \footnote{Parameters are variables 
+that a function can use but are defined elsewhere. The real values of the variables are passed on to the funciton in the call.} in \autoref{alg:stream1v1} are the 
+things that change during the execution of the model or that are calculated beforehand and not constants. $S$ for instance is not constant (well at this point it is but in \autoref{sec:daynight} 
+we change that) amd the current time is obviously not constant. All constants can be found in \autoref{sec:cp}.
 
 \begin{algorithm}[hbt]
     \SetAlgoLined
-    $\delta t \leftarrow 60 \cdot 5$ \;
+    \SetKwInput{Input}{Input}
-    $\sigma \leftarrow 5.67 \cdot 10^{-8}$ \;
+    \SetKwInOut{Output}{Output}
-    $\epsilon \leftarrow 0.75$ \;
+    \Input{time $t$, amount of energy that hits the planet $S$}
-    $C_p \leftarrow 10^7$ \;
+    \Output{Temperature of the planet $T_p$, temperature of the atmosphere $T_a$}
-    $C_a \leftarrow 10^6$ \;
+    \For{$lat \in [-90, 90]$}{
-    $S \leftarrow 1370$ \;
+        \For{$lon \in [0, 360]$}{
-    $R \leftarrow 6.4 \cdot 10^6$ \;
+            $T_p[lat, lon] \leftarrow T_p[lat, lon] + \frac{\delta t (S + 4\epsilon \sigma (T_a[lat, lon])^4 - 4\sigma (T_p[lat, lon])^4)}{C_p}$ \;
-    $t \leftarrow 0$ \;
+            $T_a[lat, lon] \leftarrow T_a[lat, lon] + \frac{\delta t (\sigma (T_p[lat, lon])^4 - 2\epsilon\sigma (T_a[lat, lon])^4)}{C_a}$ \;
-
-    \While{\texttt{TRUE}}{
-        \For{$lat \in [-90, 90]$}{
-            \For{$lon \in [0, 360]$}{
-                $T_p[lat, lon] \leftarrow T_p[lat, lon] + \frac{\delta t (S + 4\epsilon \sigma (T_a[lat, lon])^4 - 4\sigma (T_p[lat, lon])^4)}{C_p}$ \;
-                $T_a[lat, lon] \leftarrow T_a[lat, lon] + \frac{\delta t (\sigma (T_p[lat, lon])^4 - 2\epsilon\sigma (T_a[lat, lon])^4)}{C_a}$ \;
-                $t \leftarrow t + \delta t$ \;
-            }
         }
     }
-    \caption{The main loop of the temperature calculations}
+    \caption{The main function of the temperature calculations}
     \label{alg:stream1v1}
 \end{algorithm}
 
-\subsection{Day/Night Cycle}
+\subsection{Day/Night Cycle} \label{sec:daynight}
-As you can see, the amount of energy that reaches the atmopsphere is constant. However this varies based on the position of the sun relative to the planet. To fix this, we have to assign a function 
+As you can see, the amount of energy that reaches the atmopsphere is constant. However this varies based on the position of the sun relative to the planet. To fix this, we have to assign a 
-to $S$ that gives the correct amount of energy that lands on that part of the planet surface. This is done in \autoref{alg:solar}. In this algorithm the term insolation is mentioned, which is $S$
+function to $S$ that gives the correct amount of energy that lands on that part of the planet surface. This is done in \autoref{alg:solar}. In this algorithm the term insolation is mentioned, 
-used in the previous formulae if you recall. We use the $\cos$ function here to map the strength of the sun to a number between $0$ and $1$. The strength is dependent on the latitude, but since 
+which is $S$ used in the previous formulae if you recall. We use the $\cos$ function here to map the strength of the sun to a number between $0$ and $1$. The strength is dependent on the latitude, 
-that is in degrees and we need it in radians we transform it to radians by multiplying it by $\frac{\pi}{180}$. This function assumes the sun is at the equinox (center of the sun is directly 
+but since that is in degrees and we need it in radians we transform it to radians by multiplying it by $\frac{\pi}{180}$. This function assumes the sun is at the equinox (center of the sun is 
-above the equator) \cite{equinox} at at all times. The second $\cos$ is needed to simulate the longitude that the sun has moved over the longitude of the equator. For that we need the difference 
+directly above the equator) \cite{equinox} at at all times. The second $\cos$ is needed to simulate the longitude that the sun has moved over the longitude of the equator. For that we need the 
-between the longitude of the point we want to calculate the energy for, and the longitude of the sun. The longitude of the sun is of course linked to the current time (as the sun is in a different\
+difference between the longitude of the point we want to calculate the energy for, and the longitude of the sun. The longitude of the sun is of course linked to the current time (as the sun is 
-position at 5:00 than at 15:00). So we need to map the current time in seconds to the interval $[0,$ seconds in a day$]$. Therefore we need the mod function. The mod function works like this:
+in a different position at 5:00 than at 15:00). So we need to map the current time in seconds to the interval $[0,$ seconds in a day$]$. Therefore we need the mod function. The mod function 
-$x$ mod $y$ means subtract all multiples of $y$ from $x$ such that $0 \leq x < y$. So to map the current time to a time within one day, we do $t$ mod $d$ where $t$ is the current time and $d$ is 
+works like this: $x$ mod $y$ means subtract all multiples of $y$ from $x$ such that $0 \leq x < y$. So to map the current time to a time within one day, we do $-t$ mod $d$ where $-t$ is the 
-the amount of seconds in a day. When we did the calculation specified in \autoref{alg:solar} we return the final value (which means that the function call is "replaced" \footnote{Replaced is not 
+current time and $d$ is the amount of seconds in a day. We need $-t$ as this ensures that the sun moves in the right direction, with $t$ the sun would move in the opposite direction in our model 
-necessarily the right word, it is more like a mathematical function $f(x)$ where $y = f(x)$. You give it an $x$ and the value that correpsonds to that $x$ is saved in $y$. So you can view the 
+than how it would move in real life. When we did the calculation specified in \autoref{alg:solar} we return the final value (which means that the function call is "replaced" 
-function call in pseudocode as a value that is calculated by a different function which is then used like a regular number.} by the value that the function calculates). If the final value is less 
+\footnote{Replaced is not necessarily the right word, it is more like a mathematical function $f(x)$ where $y = f(x)$. You give it an $x$ and the value that correpsonds to that $x$ is saved in 
-than 0, we need to return 0 as the sun cannot suck energy out of the planet (that it does not radiate itself, which would happen if a negative value is returned).
+$y$. So you can view the function call in pseudocode as a value that is calculated by a different function which is then used like a regular number.} by the value that the function calculates). 
-
+If the final value is less than 0, we need to return 0 as the sun cannot suck energy out of the planet (that it does not radiate itself, which would happen if a negative value is returned).
-In the second stream, it was revealed that $t$ mod $d$ in \autoref{alg:solar} should be $-t$ mod $d$ such that the sun moves in the right direction. In the first stream the sun would move to the 
-right (west to east), however the sun moves to the left (east to west) and so the time must be flipped in order for the model to be correct.
 
 \begin{algorithm}[hbt]
     \SetAlgoLined
@@ -187,42 +185,20 @@ right (west to east), however the sun moves to the left (east to west) and so th
 \end{algorithm}
 
 By implementing \autoref{alg:solar}, \autoref{alg:stream1v1} must be changed as well, as $S$ is no longer constant for the whole planet surface. So let us do that in \autoref{alg:stream1v2}. Note 
-that $S$ is defined as the call to \autoref{alg:solar} (as is showcased by the text \texttt{solar}). In case you are unfamiliar with calls, defining a function is defining how it works and 
+that $S$ is defined as the call to \autoref{alg:solar}. 
-calling a function is actually using it.
-
 \begin{algorithm}[hbt]
     \SetAlgoLined
-    $\delta t \leftarrow 60 \cdot 5$ \;
+    \SetKwInput{Input}{Input}
-    $\sigma \leftarrow 5.67 \cdot 10^{-8}$ \;
+    \SetKwInOut{Output}{Output}
-    $\epsilon \leftarrow 0.75$ \;
+    \Input{time $t$, amount of energy that hits the planet $S$}
-    $C_p \leftarrow 10^7$ \;
+    \Output{Temperature of the planet $T_p$, temperature of the atmosphere $T_a$}
-    $C_a \leftarrow 10^7$ \;
-    $I \leftarrow 1370$ \;
-    $R \leftarrow 6.4 \cdot 10^6$ \;
-    $t \leftarrow 0$ \;
-    $day \leftarrow 60 \cdot 60 \cdot 24$ \;
-    $S \leftarrow$ \texttt{solar($I$, $lat$, $lon$, $t$, $day$)} \;
-    $nlat$ is the amount of latitude points in the interval $[0, 90]$, how you divide them is your own choice. \;
-    $nlot$ is the amount of longitude points in the interval $[0, 360]$, how you divide them is your own choice. \;
 
-    \While{\texttt{TRUE}}{
+    \For{$lat \in [-nlat, nlat]$}{
-        \For{$lat \in [-nlat, nlat]$}{
+        \For{$lon \in [0, nlot]$}{
-            \For{$lon \in [0, nlot]$}{
+            $T_p[lat, lon] \leftarrow T_p[lat, lon] + \frac{\delta t (S + 4\epsilon \sigma (T_a[lat, lon])^4 - 4\sigma (T_p[lat, lon])^4)}{C_p}$ \;
-                $T_p[lat, lon] \leftarrow T_p[lat, lon] + \frac{\delta t (S + 4\epsilon \sigma (T_a[lat, lon])^4 - 4\sigma (T_p[lat, lon])^4)}{C_p}$ \;
+            $T_a[lat, lon] \leftarrow T_a[lat, lon] + \frac{\delta t (\sigma (T_p[lat, lon])^4 - 2\epsilon\sigma (T_a[lat, lon])^4)}{C_a}$ \;
-                $T_a[lat, lon] \leftarrow T_a[lat, lon] + \frac{\delta t (\sigma (T_p[lat, lon])^4 - 2\epsilon\sigma (T_a[lat, lon])^4)}{C_a}$ \;
-                $t \leftarrow t + \delta t$ \;
-            }
         }
     }
-    \caption{The main loop of the temperature calculations}
+    \caption{The main function of the temperature calculations}
     \label{alg:stream1v2}
-\end{algorithm}
+\end{algorithm}
-
-\autoref{alg:stream1v2} calculates the values that are plotted (which is not discussed here as that is Python specific). Due to the \texttt{WHILE(TRUE)} loop, this calculation never finishes and
-allows us to simulate days, weeks, months and even years of heat exchange all conveniently plotted in a graph. In Simon's implementation, the graphs update in realtime, meaning that whenever a 
-round of calculations has finished, they are immediately processed to be displayed in the graph. 
-
-However other forms of looking at the calculated data can be implemented, like writing a table to a txt file, saving the generated grpahs at a certain interval or spewing all the data into a csv 
-dataset. The possibilities are endless, and the whole goal of the model is for it to be modular. Meaning that if you want to do something with it (like have a multi-layered atmosphere instead of 
-a single layer atmosphere) you can just write some lines of code and run the model and it should still work. Therefore you can write your own extensions of the model to fit it to your needs and 
-requirements.

tex-docs/topics/util_funcs.tex

@@ -0,0 +1,192 @@
+\section{Utility Functions}
+With the control panel defined and explained, let us move over to some utility functions. Functions that can be used in all kinds of calculations, which we might need more often. In general it 
+concerns functions like calculating the gradient, the lacplacian or interpolation. 
+
+\subsection{Gradients}
+Let us define the gradient in the $x, y$ and $z$ directions. The functions can be found in \autoref{alg:gradient x}, \autoref{alg:gradient y} and \autoref{alg:gradient z}. We use these functions 
+in various other algorithms as the gradient (also known as derivative) is often used in physics. It denotes the rate of change, how much something changes over time. Velocity for instance denotes
+how far you move in a given time. Which is a rate of change, how much your distance to a given point changes over time. 
+
+In \autoref{alg:gradient z} $a.dimensions$ is the attribute that tells us how deeply nested the array $a$ is. If the result is $1$ we have just a normal array, if it is $2$ we have a double array 
+(an array at each index of the array) which is also called a matrix and if it is $3$ we have a triple array. We need this because we have a one-dimensional case, for when we do not use multiple 
+layers and a three-dimensional case for when we do use multiple layers. This distinction is needed to avoid errors being thrown when running the model with one or multiple layers. 
+
+This same concept can be seen in \autoref{alg:gradient x} and \autoref{alg:gradient y}, though here we check if $k$ is defined or \texttt{NULL}. We do this as sometimes we want to use this 
+function for matrices that does not have the third dimension. Hence we define a default value for $k$ which is \texttt{NULL}. \texttt{NULL} is a special value in computer science. It represents 
+nothing. This can be useful sometimes if you declare a variable to be something but it is referring to something that has been deleted or it is returned when some function fails. It usually 
+indicates that something special is going on. So here we use it in the special case where we do not want to consider the third dimension in the gradient. We also use forward differencing 
+(calculating the gradient by taking the difference of the cell and the next/previous cell, multiplied by $2$ to keep it fair) in \autoref{alg:gradient y} as that gives better results for the 
+calculations we will do later on.
+
+\begin{algorithm}[hbt]
+    \SetKwInOut{Input}{Input}
+    \SetKwInOut{Output}{Output}
+    \Input{Matrix (double array) $a$, first index $i$, second index $j$, third index $k$ with default value \texttt{NULL}}
+    \Output{Gradient in the $x$ direction}
+    \eIf{$k == \texttt{NULL}$}{
+        $grad \leftarrow \frac{a[i, (j + 1)\text{ mod } nlon] - a[i, (j - 1) \text{ mod } nlon]}{\delta x[i]}$ \;
+    }{
+        $grad \leftarrow \frac{a[i, (j + 1)\text{ mod } nlon, k] - a[i, (j - 1) \text{ mod } nlon, k]}{\delta x[i]}$ \;
+    }
+    \Return{$grad$} \;
+    \caption{Calculating the gradient in the $x$ direction}
+    \label{alg:gradient x}
+\end{algorithm}
+
+\begin{algorithm}[hbt]
+    \SetKwInOut{Input}{Input}
+    \SetKwInOut{Output}{Output}
+    \Input{Matrix (double array) $a$, first index $i$, second index $j$, third index $k$ with default value \texttt{NULL}}
+    \Output{Gradient in the $y$ direction}
+    \eIf{$k == \texttt{NULL}$}{
+        \uIf{$i == 0$}{
+            $grad \leftarrow 2 \frac{a[i + 1, j] - a[i, j]}{\delta y}$ \;
+        }\uElseIf{$i == nlat - 1$}{
+            $grad \leftarrow 2 \frac{a[i, j] - a[i - 1, j]}{\delta y}$ \;
+        }\uElse{
+            $grad \leftarrow \frac{a[i + 1, j] - a[i - 1 j]}{\delta y}$ \;
+        }
+    }{
+        \uIf{$i == 0$}{
+            $grad \leftarrow 2 \frac{a[i + 1, j, k] - a[i, j, k]}{\delta y}$ \;
+        }\uElseIf{$i == nlat - 1$}{
+            $grad \leftarrow 2 \frac{a[i, j, k] - a[i - 1, j, k]}{\delta y}$ \;
+        }\uElse{
+            $grad \leftarrow \frac{a[i + 1, j] - a[i - 1 j]}{\delta y}$ \;
+        }
+    }
+    \Return $grad$ \;
+    \caption{Calculating the gradient in the $y$ direction}
+    \label{alg:gradient y}
+\end{algorithm}
+
+\begin{algorithm}[hbt]
+    \SetKwInOut{Input}{Input}
+    \SetKwInOut{Output}{Output}
+    \Input{Matrix (double array) $a$, first index $i$, second index $j$, third index $k$}
+    \Output{Gradient in the $z$ direction}
+    \uIf{$a.dimensions == 1$}{
+        \uIf{$k == 0$}{
+            $grad \leftarrow \frac{a[k + 1] - a[k]}{\delta z[k]}$ \;
+        }\uElseIf{$k == nlevels - 1$}{
+            $grad \leftarrow \frac{a[k] - a[k - 1]}{\delta z[k]}$ \;
+        }\uElse{
+            $grad \leftarrow \frac{a[k + 1] - a[k - 1]}{2\delta z[k]}$ \;
+        }
+    } \uElse {
+        \uIf{$k == 0$}{
+            $grad \leftarrow \frac{a[i, j, k + 1] - a[i, j, k]}{\delta z[k]}$ \;
+        }\uElseIf{$k == nlevels - 1$}{
+            $grad \leftarrow \frac{a[i, j, k] - a[i, j, k - 1]}{\delta z[k]}$ \;
+        }\uElse{
+            $grad \leftarrow \frac{a[i, j, k + 1] - a[i, j, k - 1]}{2\delta z[k]}$ \;
+        }
+    }
+    
+    \Return $grad$ \;
+    \caption{Calculating the gradient in the $z$ direction}
+    \label{alg:gradient z}
+\end{algorithm}
+
+\subsection{Laplacian Operator} \label{sec:laplace}
+The Laplacian operator ($\nabla^2$, sometimes also seen as $\Delta$) has two definitions, one for a vector field and one for a scalar field. The two concepts are not indpendent, a vector field 
+is composed of scalar fields \cite{vectorscalarfields}. Let us define a vector field first. A vector field is a function whose domain and range are a subset of the Eucledian $\mathbb{R}^3$ space. 
+A scalar field is then a function consisting out of several real variables (meaning that the variables can only take real numbers as valid values). So for instance the circle equation 
+$x^2 + y^2 = r^2$ is a scalar field as $x, y$ and $r$ are only allowed to take real numbers as their values. 
+
+With the vector and scalar fields defined, let us take a look at the Laplacian operator. For a scalar field $\phi$ the laplacian operator is defined as the divergence of the gradient of $\phi$
+\cite{laplacian}. But what are the divergence and gradient? The gradient is defined in \autoref{eq:gradient} and the divergence is defined in \autoref{eq:divergence}. Here $\phi$ is a vector 
+with components $x, y, z$ and $\Phi$ is a vector field with components $x, y, z$. $\Phi_1, \Phi_2$ and $\Phi_3$ refer to the functions that result in the corresponding $x, y$ and $z$ values 
+\cite{vectorscalarfields}. Also, $i, j$ and $k$ are the basis vectors of $\mathbb{R^3}$, and the multiplication of each term with their basis vector results in $\Phi_1, \Phi_2$ and $\Phi_3$
+respectively. If we then combine the two we get the Laplacian operator, as in \autoref{eq:laplacian scalar}.
+
+\begin{subequations}
+    \begin{equation}
+        \text{grad } \phi = \nabla \phi = \frac{\delta \phi}{\delta x}i + \frac{\delta \phi}{\delta y}j + \frac{\delta \phi}{\delta z}k
+        \label{eq:gradient}
+    \end{equation}
+    \begin{equation}
+        \text{div} \Phi = \nabla \cdot \Phi = \frac{\delta \Phi_1}{\delta x} + \frac{\delta \Phi_2}{\delta y} + \frac{\delta \Phi_3}{\delta z}
+        \label{eq:divergence}
+    \end{equation}
+    \begin{equation}
+        \nabla^2 \phi = \nabla \cdot \nabla \phi = \frac{\delta^2 \phi}{\delta x^2} + \frac{\delta^2 \phi}{\delta y^2} + \frac{\delta^2 \phi}{\delta z^2}
+        \label{eq:laplacian scalar}
+    \end{equation}
+\end{subequations}
+
+For a vector field $\Phi$ the Laplacian operator is defined as in \autoref{eq:laplacian vector}. Which essential boils down to taking the Laplacian operator of each function and multiply it by
+the basis vector.
+
+\begin{equation}
+    \nabla^2 \Phi = (\nabla^2 \Phi_1)i + (\nabla^2 \Phi_2)j + (\nabla^2 \Phi_3)k
+    \label{eq:laplacian vector}
+\end{equation}
+
+The new code can be found in \autoref{alg:laplacian}. $\Delta_x$ and $\Delta_y$ in \autoref{alg:laplacian} represents the calls to \autoref{alg:gradient x} and \autoref{alg:gradient y} 
+respectively.
+
+\begin{algorithm}[hbt]
+    \SetKwInOut{Input}{Input}
+    \SetKwInOut{Output}{Output}
+    \Input{A matrix (double array) a}
+    \Output{A matrix (double array) with results for the laplacian operator for each element}
+    \eIf{$a.dimensions == 2$}{
+        \For{$lat \in [1, nlat - 1]$}{
+            \For{$lon \in [0, nlon]$}{
+                $output[lat, lon] \leftarrow \frac{\Delta_x(a, lat, (lon + 1) \text{ mod } nlon) - \Delta_x(a, lat, (lon - 1) \text{ mod } nlon)}{\delta x[lat]} + \frac{\Delta_y(a, lat + 1, lon) - 
+                \Delta_y(a, lat - 1, lon)}{\delta y}$\;
+            }
+        }
+    }{
+        \For{$lat \in [1, nlat - 1]$}{
+            \For{$lon \in [0, nlon]$}{
+                \For{$k \in [0, nlevels - 1]$}{
+                    $output[lat, lon, k] \leftarrow \frac{\Delta_x(a, lat, (lon + 1) \text{ mod } nlon, k) - \Delta_x(a, lat, (lon - 1) \text{ mod } nlon, k)}{\delta x[lat]} + \frac{\Delta_y(a, 
+                    lat + 1, lon, k) - \Delta_y(a, lat - 1, lon, k)}{\delta y} + \frac{\Delta_z(a, lat, lon, k + 1) - \Delta_z(a, lat, lon, k + 1)}{2\delta z[k]}$\;
+                }
+            }
+        }
+    }
+    
+    \Return{$ouput$} \;
+    \caption{Calculate the laplacian operator over a matrix a}
+    \label{alg:laplacian}
+\end{algorithm}
+
+\subsection{Divergence}
+As we expect to use the divergence operator more often throughout our model, let us define a seperate function for it in \autoref{alg:divergence}. $\Delta_x$ and $\Delta_y$ in 
+\autoref{alg:divergence} represents the calls to \autoref{alg:gradient x} and \autoref{alg:gradient y} respectively. We do the multiplication with the velocity vectors $u, v$ and $w$ here already, 
+as we expect that we might use it in combination with the divergence operator more frequently. What those vectors are and represent we will discuss in %insert velocity reference here
+
+\begin{algorithm}[!hbt]
+    \SetKwInOut{Input}{Input}
+    \SetKwInOut{Output}{Output}
+    \Input{A matrix (double array) $a$}
+    \Output{A matrix (double array) containing the result of the divergence operator taken over that element}
+    $dim_1 \leftarrow \text{ Length of } a \text{ in the first dimension}$ \;
+    \For{$i \in [0, dim_1]$}{
+        $dim_2 \leftarrow \text{ Length of } a \text{ in the second dimension (i.e. the length of the array stored at index } i)$ \;
+        \For{$j \in [0, dim_2]$}{
+            $dim_3 \leftarrow \text{ Length of } a \text{ in the third dimension}$ \;
+            \For{$k \in [0, dim_3]$}{
+                $output[i, j, k] \leftarrow \Delta_x(au, i, j, k) + \Delta_y(av, i, j, k) + \Delta_z(aw, i, j, k)$ \;
+            }
+        }
+    }
+    \Return{$output$} \;
+    \caption{Calculate the result of the divergence operator on a vector}
+    \label{alg:divergence}
+\end{algorithm}
+
+\subsection{Interpolation} \label{sec:interpolation}
+Interpolation is a form of estimation, where one has a set of data points and desires to know the values of other data points that are not in the original set of data points\cite{interpolation}. 
+Based on the original data points, it is estimated what the values of the new data points will be. There are various forms of interpolation like linear interpolation, polynomial interpolation 
+and spline interpolation. The CLAuDE model uses linear interpolation which is specified in \autoref{eq:interpolation}. Here $z$ is the point inbetween the known data points $x$ and $y$. 
+$\lambda$ is the factor that tells us how close $z$ is to $y$ in the interval $[0, 1]$. If $z$ is very close to $y$, $\lambda$ will have the value on the larger end of the interval, like 0.9.
+Whereas if $z$ is close to $x$ then $\lambda$ will have a value on the lower end of the interval, like 0.1.
+
+\begin{equation}
+    z = (1 - \lambda)x + \lambda y
+    \label{eq:interpolation}
+\end{equation}