\section{The Beginning}
\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 
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 $\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 
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}
    \item $\Delta U$: The change of internal energy ($J$) \cite{thermo1}.
    \item $C$: The specific heat capacity of the object, i.e. how much energy is required to heat the object by one degree Kelvin ($\frac{J}{K}$).
    \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 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.

\begin{subequations}
    \begin{equation}
        S = SB = \sigma T^4
        \label{eq:stefan-boltzmann}
    \end{equation}
    \begin{equation}
        S - \sigma T^4 = \Delta U = C \Delta T
        \label{eq:sb rewritten}
    \end{equation}
    \begin{equation}
        \Delta T = \frac{\delta t(S - \sigma T^4)}{C}
        \label{eq:sb rewritten2}
    \end{equation}
    \label{eq:basis}
\end{subequations}

The set of equations in \autoref{eq:basis} form the basis of the temperature exchange of the planet. However two crucial aspects are missing. Only half of the planet will be receiving light from 
the sun at once, and the planet is a sphere. So we need to account for both in our equation. We do that in \autoref{eq:basis sphere correction}. We view the energy reacing the atmosphere as a 
circular area of energy, with the equation for the are of a circle being \autoref{eq:basis circle} \cite{areaCircle}. The area of a sphere is in \autoref{eq:basis sphere} \cite{areaSphere}. In 
both equations, $r$ is the radius of the circle/sphere. By using \autoref{eq:basis circle} and \autoref{eq:basis sphere} in \autoref{eq:sb rewritten2} we get \autoref{eq:basis sphere2} where 
$r$ is replaced by $R$. It is common in physics literature to use capitals for large objects like planets. However we are not done yet since we can divide some stuff out. We end up with 
\autoref{eq:basis sphere final} as the final equation we are going to use.

\begin{subequations}
    \begin{equation}
        \pi r^2
        \label{eq:basis circle}
    \end{equation}
    \begin{equation}
        4 \pi r^2
        \label{eq:basis sphere}
    \end{equation}
    \begin{equation}
        \Delta T = \frac{\delta t (\pi R^2S - 4\pi R^2\sigma T^4)}{4\pi CR^2}
        \label{eq:basis sphere2}
    \end{equation}
    \begin{equation}
        \Delta T = \frac{\delta t (S - 4\sigma T^4)}{4C}
        \label{eq:basis sphere final}
    \end{equation}
    \label{eq:basis sphere correction}
\end{subequations}

\subsection{Insolation}

With the current equation we calculate the global average surface temperature of the planet itself. However, this planet does not have an atmosphere just yet. Basically we modelled the 
temperature of a rock floating in space, let's change that with \autoref{eq:atmos}. Here we assume that the area of the atmosphere is equal to the area of the planet surface. Obviously
this assumption is false, as the atmosphere is a sphere that is larger in radius than the planet, however the difference is not significant enough to account for it. We also define the
atmosphere as a single layer. This is due to the accessibility of the model, we want to make it accessible, not university simulation grade. One thing to take into account for the 
atmosphere is that it only partially absorbs energy. The sun (or a similar star) is relatively hot and sends out energy waves (radiation) with relatively low wavelengths. The planet is 
relatively cold and sends out energy at long wavelengths. As a side note, all objects radiate energy. You can verify this by leaving something in the sun on a hot day for a while and 
almost touch it later. You can feel the heat radiating from the object. The planet is no exception and radiates heat as well, though at a different wavelength than the sun. The 
atmosphere absorbs longer wavelengths better than short wavelengths. For simplicity's sake we say that all of the sun's energy does not get absorbed by the atmosphere. The planet's 
radiation will be absorbed partially by the atmosphere. Some of the energy that the atmosphere absorbs is radiated into space and some of that energy is radiated back onto the planet's 
surface. We need to adjust \autoref{eq:basis sphere final} to account for the energy being radiated from the atmosphere back at the planet surface.

So let us denote the temperature of the planet surface as $T_p$ and the temperature of the atmosphere as $T_a$. Let us also write the specific heat capacity of the planet surface as $C_p$ 
instead of $C$. We add the term in \autoref{eq:atmos on surface improved} to \autoref{eq:basis sphere final} in \autoref{eq:surface change}. In \autoref{eq:atmos on surface}, $\epsilon$ is the 
absorbtivity of the atmosphere, the fraction of energy that the atmosphere absorbs. We divided \autoref{eq:atmos on surface} by $\pi R^2$ as we did that with \autoref{eq:basis sphere final} as 
well, so we needed to make it match that division.

\begin{subequations}
    \begin{equation}
        4\pi R^2 \epsilon \sigma T_a^4
        \label{eq:atmos on surface}
    \end{equation}
    \begin{equation}
        4\epsilon \sigma T_a^4
        \label{eq:atmos on surface improved}
    \end{equation}
    \begin{equation}
        \Delta T_p = \frac{\delta t (S + 4\epsilon \sigma T_a^4 - 4\sigma T_p^4)}{4C_p}
        \label{eq:surface change}
    \end{equation}
    \label{eq:atmos}
\end{subequations}

As you probably expected, the atmosphere can change in temperature as well. This is modelled by \autoref{eq:atmos change}, which is very similar to \autoref{eq:surface change}. There are
some key differences though. Instead of subtracting the radiated heat of the atmosphere once we do it twice. This is because the atmosphere radiates heat into space and towards the 
surface of the planet, which are two outgoing streams of energy instead of one for the planet (as the planet obviously cannot radiate energy anywhere else than into the atmosphere).
$C_a$ is the specific heat capacity of the atmosphere.

\begin{equation}
    \Delta T_a = \frac{\delta t (\sigma T_p^4 - 2\epsilon\sigma T_a^4)}{C_a}
    \label{eq:atmos change}
\end{equation}

\subsection{The Latitude Longitude Grid}
With the current model, we only calculate the global average temperature. To calculate the temperature change along the surface and atmosphere at different points, we are going to use a grid.
Fortunately the world has already defined such a grid for us, the latitude longitude grid \cite{latlong}. The latitude is the coordinate running from the south pole to the north pole, with -90 
being the south pole and 90 being the north pole. The longitude runs parallel to the equator and runs from 0 to 360 which is the amount of degrees that an angle can take when calculating the 
angle of a circle. So 0 degrees longitude is the same place as 360 degrees longitude. To do this however we need to move on from mathematical formulae to code (or in this case pseudocode).

Pseudocode is a representation of real code. It is meant to be an abstraction of code such that it does not matter how you present it, but every coder should be able to read it and implement 
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 
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 
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.
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 
individual cell of that grid indexed by a specific latitude longitude combination.

\begin{algorithm}[hbt]
    \SetAlgoLined
    $\delta t \leftarrow 60 \cdot 5$ \;
    $\sigma \leftarrow 5.67 \cdot 10^{-8}$ \;
    $\epsilon \leftarrow 0.75$ \;
    $C_p \leftarrow 10^7$ \;
    $C_a \leftarrow 10^6$ \;
    $S \leftarrow 1370$ \;
    $R \leftarrow 6.4 \cdot 10^6$ \;
    $t \leftarrow 0$ \;

    \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}
    \label{alg:stream1v1}
\end{algorithm}

\subsection{Day/Night Cycle}
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 
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$
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 
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 
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 
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\
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:
$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 
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 
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 
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
    \SetKwInput{Input}{Input}
    \SetKwInOut{Output}{Output}
    \Input{insolation $ins$, latitude $lat$, longitude $lon$, time $t$, time in a day $d$}
    \Output{Amount of energy $S$ that hits the planet surface at the given latitude-time combination.}
    $longitude \leftarrow 360 \cdot \frac{(-t \text{ mod } d)}{d}$ \;
    $S \leftarrow ins \cdot \cos(lat \frac{\pi}{180}) \cos((lon - longitude) \cdot \frac{\pi}{180})$ \;
    \eIf{$S < 0$}{
        \Return{$0$}
    }{
        \Return{$S$}
    }
    \caption{Calculating the energy from the sun (or similar star) that reaches a part of the planet surface at a given latitude and time}
    \label{alg:solar}
\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 
calling a function is actually using it.

\begin{algorithm}[hbt]
    \SetAlgoLined
    $\delta t \leftarrow 60 \cdot 5$ \;
    $\sigma \leftarrow 5.67 \cdot 10^{-8}$ \;
    $\epsilon \leftarrow 0.75$ \;
    $C_p \leftarrow 10^7$ \;
    $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{$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_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}
    \label{alg:stream1v2}
\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.