Added the latest stream of 26-08-2020, where we got a working model! Hooray!

tex-docs/CLAuDE.tex

@@ -71,6 +71,8 @@ knowing).
 
 \input{streams/Stream9.tex}
 
+\input{streams/Stream10.tex}
+
 \newpage
 \input{streams/TTNMETAF.tex}
 

tex-docs/streams/Stream10.tex

@@ -0,0 +1,81 @@
+\section{Putting Our Homemade Climate Model Through Its Paces}
+Big stream this stream as we got some working code! Always great when stuff works. This stream, we tackled the radiation problem, added axial tilt to the planet, fixed vertical motion (but not
+advection), added stratospheric heating and some other code clean up stuff. This means that the big rework is getting closer and closer. How exciting!
+%TESTSSSS WE GOT WORKING CODE! HELLO WORLD!
+
+\subsection{Fixing Up the Code}
+First thing to mention is that vertical advection is still broken. Why? Because the gradient in the $z$ direction is broken. This is due to finite differencing on an exponential function. The way
+we calculate the differenc from one layer to the other is by differencing them (subtracting) which is always finite. Therefore we always get some inaccuracies. Usually that is fine, but with an 
+exponential function the differences, you guessed it, become exponentially wrong. As such, the function would eventually be so far off that the model would blow up. So we need to fix it. To 
+prevent a blow up, we have disabled the call to the gradient $z$ funciton in \autoref{alg:divergence layer}. This ensures that the horizontal bits still work, but the vertical stuff does not.
+As always, we will try to fix this in a future stream. 
+
+We also fixed up the radiation scheme, as shown in \autoref{alg:optical depth}. Basically we had the definition of $U[k + 1] = \text{something something} U[k + 1]$. This means that the definition
+was relying on itself, which is obviously impossible and wrong. So we changed it to it's current form and it is fixed, hooray!
+
+Vertical motion has also been fixed, as shown in \autoref{alg:velocity}. Due to some error in the representation of the vertical motion it did not work. So we changed from that representation to 
+another. Now the vertical velocity is proportional to the rate of change of the pressure which does work like it should.
+
+\subsection{Tilting the Planet}
+In order to model a planet that has seasons, like Earth, we need to tilt the planet. This has as effect that the sun is not always directly above the equator but is above a certain band around
+the equator as the year moves on. This means that some hemispheres receive more/less sun based on what part of the year it is. Which corresponds to the various seasons we have on Earth. But in
+order to do that, we have to change the \texttt{solar} function. The new version as shown in \autoref{alg:solar tilt} will replace \autoref{alg:solar}. Here $\alpha$ is the tilt in degrees.
+
+\begin{algorithm}
+    \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, longitude and time combination.}
+    $sun\_lon \leftarrow -t \text{ mod } d$ \;
+    $sun\_lon \leftarrow sun\_lon \cdot \frac{360}{d}$ \;
+    $sun\_lat \leftarrow \alpha\cos(\frac{2t\pi}{year})$ \;
+    $S \leftarrow insolation\cos(\frac{\pi(lat - sun\_lat)}{180})$ \;
+
+    \uIf{$S < 0$}{
+        \Return{$0$} \;
+    } \uElse {
+        $lon\_diff \leftarrow lon - sun\_lon$ \;
+        $S \leftarrow S\cos(\frac{lon\_diff\pi}{180})$ \;
+
+        \uIf{$S < 0$}{
+            \uIf{$lat + sun\_lat > 90$ or $lat + sun\_lat < -90$}{
+                \Return{$insolation\cos(\frac{\pi(lat + sun\_lat)}{180})\cos(\frac{lon\_diff\pi}{180})$} \;
+            } \uElse {
+                \Return{$0$} \;
+            }
+        } \uElse {
+            \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 tilt}
+\end{algorithm}
+
+What the code in \autoref{alg:solar tilt} does boils down to calculating the latitude and longitude of the sun and checking whether the planet receives any energy. If not return $0$ immediately.
+If so we check if the difference between the sun's longitude and the planet's longitude and calculate how much energy would hit the planet given that the sun is not straight above the equator. 
+We do this by multiplying the energy it would receive from the sun if it were above the equator $S$ by the cosine of the difference in longitudes, which represents the tilt. Then we check again 
+if the planet is receiving energy, if not we check if it happens around the poles. We do this because due to the tilt it can be the case that at certain points in the year the pole is in constant
+sunlight, i.e. the sun does not go down. This creates a sort of overshoot which needs to be accounted for. If it does this then we add the latitudes of the sun and the planet together and use
+that to calculate the energy that would hit that spot. If it is not the case that we are around the poles and we do not receive energy, then we return $0$. If it happens to be that we do receive 
+energy (so no negative values) then we return $S$.
+
+\subsection{Adding In Some Ozone (Or Something Else That Approximates It)}
+Adding in ozone in the stratosphere is hella complicated, so we leave that as an exercise to the reader as in true academic fashion. Just joking, if you want you can work on implementing ozone 
+however we opt not to because it is quite complicated. Instead we approximate it, which is decent enough for our purpose. We need to do it in \autoref{alg:optical depth} as we need to adjust the
+$Q$. We add in a check to see if we are currently calculating the radiation in the stratosphere. If so we add some radiation extra to replicate the effect of ozone. As can be seen in 
+\autoref{alg:ozone}, where we only focus on the $Q$ part of \autoref{alg:optical depth}, we add in some extra radiation based on how high the current layer calculation is, which scales with the
+height. 
+
+\begin{algorithm}
+    \For{$level \in [0, nlevels]$}{
+        $Q[level] \leftarrow - \frac{S_z(U - D, 0, 0, level)}{10^3 \cdot densityProfile[level]}$ \;
+        \uIf{$heights[level] > 20 \cdot 10^3$}{
+            $Q[level] \leftarrow Q[level] + \texttt{solar}(5, lat, lon, t) \frac{24 \cdot 60 \cdot 60(\frac{heights[level] - 20 \cdot 10^3}{10^3})^2}{30^2}$ \;
+        }
+    }
+    \caption{Replicating the effect of ozone}
+    \label{alg:ozone}
+\end{algorithm}
+
+It is at this point that we reached the state that CLAuDE is in a testable state. This means that we have the model working in such a way that we can do some simple experiments like altering how
+long a day is, what would happen if the sun would send out more energy (which usually means that it is bigger) or what would happen if you tidally lock a planet (stop it rotating completely). 

tex-docs/streams/Stream5.tex

@@ -229,7 +229,9 @@ We also need to initialise the $\epsilon$ value for each layer. We do that in \a
 \end{algorithm}
 
 Now we need to add vertical winds, or in other words add the $w$ component of the velocity vectors. We do that by editing \autoref{alg:stream3}. We change it to \autoref{alg:velocity}. Here we 
-use gravity ($g$) instead of the coriolis force ($f$) and calculate the pressure gradient in the $z$ direction.
+use gravity ($g$) instead of the coriolis force ($f$) and calculate the change in pressure. Therefore we need to store a copy of the pressure before we do any calculations. This needs to be a
+copy due to aliasing \footnote{Aliasing is assigning a different name to a variable, while it remains the same variable. Take for instance that we declare a variable $x$ and set it to be $4$. 
+Then we say $y \leftarrow x$, which you might think is the same as saying they $y \leftarrow 4$ but behind the screen it is pointing to $x$. So if $x$ changes, then so does $y$.}
 
 \begin{algorithm}
     $S_{xu} \leftarrow \texttt{gradient\_x}(u, lan, lon)$ \;
@@ -238,17 +240,18 @@ use gravity ($g$) instead of the coriolis force ($f$) and calculate the pressure
     $S_{yv} \leftarrow \texttt{gradient\_y}(v, lan, lon)$ \;
     $S_{px} \leftarrow \texttt{gradient\_x}(p, lan, lon)$ \;
     $S_{py} \leftarrow \texttt{gradient\_y}(p, lan, lon)$ \;
-    $S_{pz} \leftarrow \texttt{gradient\_z}(p, lan, lon)$ \;
     \While{\texttt{TRUE}}{
         \For{$lat \in [1, nlat - 1]$}{
             \For{$lon \in [0, nlon]$}{
                 \For{$layer \in [0, nlevels]$}{
                     $u[lan, lon, layer] \leftarrow u[lat, lon, layer] + \delta t \frac{-u[lat, lon, layer]S_{xu} - v[lat, lon, layer]S_{yu} + f[lat]v[lat, lon, layer] - S_{px}}{\rho}$ \;
                     $v[lan, lon, layer] \leftarrow v[lat, lon, layer] + \delta t \frac{-u[lat, lon, layer]S_{xv} - v[lat, lon, layer]S_{yv} - f[lat]u[lat, lon, layer] - S_{py}}{\rho}$ \;
-                    $w[lan, lon, layer] \leftarrow w[lat, lon, layer] + \delta t (\frac{S_{pz}}{\rho[lat, lon, layer]} + g)$ \;
+                    $w[lan, lon, layer] \leftarrow w[lat, lon, layer] - \frac{p[lat, lon, layer] - p_o[lat, lon, layer]}{\delta t\rho[lat, lon, layer]g}$ \;
                 }
             }
         }
+
+        $p_o \leftarrow copy(p)$ \;
     }
     \caption{Calculating the flow of the atmosphere (wind)}
     \label{alg:velocity}

tex-docs/streams/Stream8.tex

@@ -84,17 +84,13 @@ $S_z$ represents the call to \autoref{alg:gradient z layer}. \texttt{solar} repr
             $U[0] \leftarrow \sigma T_p[lat, lon]^4$ \;
 
             \For{$level \in [1, nlevels]$}{
-                $U[level] \leftarrow U[level - 1] + (\tau[level] - \tau[level - 1])(U[level - 1] - \sigma \cdot (mean(T_a[:, :, level]))^4)$ \;
+                $U[level] \leftarrow U[level - 1] - \frac{(\tau[level] - \tau[level - 1])(\sigma \cdot (mean(T_a[:, :, level]))^4)}{1 + (\tau[level - 1] - \tau[level])}$ \;
             }
 
             $D[nlevels - 1] \leftarrow 0$ \;
 
             \For{$level \in [nlevels - 1, 0]$}{
-                \uIf{$level == 0$}{
+                $D[level] \leftarrow D[level + 1] - \frac{(\tau[level + 1] - \tau[level])(\sigma \cdot (mean(T_a[:, :, level]))^4)}{1 + (\tau[level] - \tau[level + 1])}$ \;
-                    $D[level] \leftarrow D[level + 1] + (\tau[level + 1] - \tau[level])(\sigma \cdot (T_a[lat, lon, level])^4) - D[level + 1]$ \;
-                }\uElse{
-                    $D[level] \leftarrow D[level + 1] + (\tau[level] - \tau[level - 1])(\sigma \cdot (T_a[lat, lon, level])^4) - D[level + 1]$ \;
-                }
             }
 
             \For{$level \in [0, nlevels]$}{