Merge pull request #6 from TechWizzart/master

Updated the manual to include the stream from 19-08-2020

tex-docs/CLAuDE.tex

@@ -47,6 +47,12 @@ One important thing to note, the layout may change significantly when new sectio
 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. 
+The fork will usually be more up to date than the version on the Planet Factory repository as Simon needs to merge pull requests into the repository. However I can update the fork freely so if a 
+particular stream is missing in the version on the Planet Factory repository, check the fork/Discord whether there is a newer version. If that is not the case, you just have to be a bit more 
+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}
@@ -63,6 +69,10 @@ manual) figure out the placement of the algorithm blocks, which may or may not b
 
 \input{streams/Stream8.tex}
 
+\input{streams/Stream9.tex}
+
+\input{streams/Stream10.tex}
+
 \newpage
 \input{streams/TTNMETAF.tex}
 

tex-docs/references.bib

@@ -15,6 +15,14 @@ editor={Reed, CaptainEditor}
 author={Clark, Simon}
 }
 
+@misc{polarPlane, 
+howpublished="\url{https://www.twitch.tv/videos/715183053}", 
+author={Clark, Simon},
+journal={Twitch}, 
+year={2020}, 
+month={Aug}
+}
+
 @misc{SI, 
 title={SI Units}, 
 howpublished="\url{https://www.nist.gov/pml/weights-and-measures/metric-si/si-units}", 
@@ -307,3 +315,29 @@ author={Frierson, Dargan M. W. and Held, Isaac M. and Zurita-Gotor, Pablo},
 year={2006}, 
 pages={2548–2566}
 }
+
+@misc{claudeGit, 
+title={Planet-Factory/claude}, 
+howpublished="\url{https://github.com/Planet-Factory/claude}", 
+journal={GitHub}, 
+author={Clark, Simon}, 
+editor={Reid, Peter and Marsden, Aren},
+translator={Baggen, Sam}
+}
+
+@misc{nomGit, 
+title={TechWizzart/claude}, 
+howpublished="\url{https://github.com/TechWizzart/claude}", 
+journal={GitHub}, 
+author={Baggen, Sam}
+}
+
+@misc{pythagoras, 
+title={Pythagorean theorem}, 
+howpublished="\url{https://en.wikipedia.org/wiki/Pythagorean_theorem}", 
+journal={Wikipedia}, 
+publisher={Wikimedia Foundation}, 
+author={Boldt, Axel}, 
+year={2020}, 
+month={Aug}
+}

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/Stream6.tex

@@ -12,7 +12,8 @@ the total advection (and because it was at this time broken). %The Laplacian has
 Another thing that we found out was broken is the vertical momentum. We tried to add it, ran into problems and just set it to 0 to fix the other problems that occured. One of those problems was 
 a wrong initialisation of the density. We basically told the model that the density is the same on every layer of the atmosphere, which is obviously not true. Hence we need to adjust that. The 
 new initialisation is described in \autoref{alg:density}. Note that the $\rho[:,: i]$ notation means that for every index in the first and second dimension, only change the value for the index $i$
-in the third dimension.
+in the third dimension. As soon as the vertical momentum has been fixed it will be fixed in the correct spot. If a lot of the code changes then it will probably be in another section and I will 
+insert a reference to that section.
 
 \begin{algorithm}
     $\rho[:, :, 0] \leftarrow 1.3$ \;

tex-docs/streams/Stream8.tex

@@ -74,7 +74,7 @@ replaces the temperature calculations we have done in \autoref{alg:temperature l
 and is located near all the other constants in the code, henceforth we will refer to this section in the code as the control panel, since most if not all of the constants can be tweaked here. 
 $\tau_0$ is a function that gives the surface optical depth for a given latitude. The corresponding equation can be found in \autoref{eq:optical depth surface} \cite{simon}. Translating this 
 into code is left as an exercise to the reader. $U[0]$ is the boundary condition discussed before (being the same as \autoref{eq:stefan-boltzmann}), just as $D[nlevels]$ is the boundary condition. 
-$S_z$ represents the call to \autoref{alg:gradient z layer}.
+$S_z$ represents the call to \autoref{alg:gradient z layer}. \texttt{solar} represents the call to \autoref{alg:solar}.
 
 \begin{algorithm}
     \For{$lat \in [-nlat, nlat]$}{
@@ -84,13 +84,13 @@ $S_z$ represents the call to \autoref{alg:gradient z layer}.
             $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] \leftarrow 0$ \;
+            $D[nlevels - 1] \leftarrow 0$ \;
 
             \For{$level \in [nlevels - 1, 0]$}{
-                $D[level] \leftarrow (\tau[level] - \tau[level - 1])(\sigma \cdot (mean(T_a[:, :, level]))^4) - D[level + 1]$ \;
+                $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])}$ \;
             }
 
             \For{$level \in [0, nlevels]$}{
@@ -98,6 +98,10 @@ $S_z$ represents the call to \autoref{alg:gradient z layer}.
             }
 
             $T_a[lat, lon, :] \leftarrow T_a[lat, lon, :] + Q$ \;
+
+            $S \leftarrow \texttt{solar}(I, lat, lon, t)$ \;
+
+            $T_p[lat, lon] \leftarrow T_p[lat, lon] \frac{\delta t((1 - a[lat, lon]) S + S_z(D, 0, 0, 0) - \sigma T_p[lat, lon]^4)}{C_p[lat ,lon]}$ \;
         }
     }
     \caption{Adding in radiation}
@@ -109,5 +113,5 @@ $S_z$ represents the call to \autoref{alg:gradient z layer}.
     \label{eq:optical depth surface}
 \end{equation}
 
-Note that in this form, it did not work on stream yet. This may be due to a coding error or to a missing equation, constant or something similar. If it turns out to be a simple fix, then it will
+%Note that in this form, it did not work on stream yet. This may be due to a coding error or to a missing equation, constant or something similar. If it turns out to be a simple fix, then it will
-be fixed in this section. If a lot of other things change in order for the fix to work, then it will probably be a seperate section with a reference to that section here.
+%be fixed in this section. If a lot of other things change in order for the fix to work, then it will probably be a seperate section with a reference to that section here.

tex-docs/streams/Stream9.tex

@@ -0,0 +1,79 @@
+\section{Starting to Deal With the Poles}
+It is time to deal with the pole situation. The north and south poles that is, not the lovely people over in Poland. We run into problems because the latitude longitude grid cells become to small 
+near the poles. Therefore, the magnitudes no longer fit into one cell and overflow into other cells which makes everything kind of funky. So we need to fix that, and we do that by a planar 
+approximation. 
+
+\subsection{The Theory Behind the Planar Approximation}
+As said earlier, the grid cells on the latitude longitude grid get closer together the closer you get to the poles which poses problems. To fix this, we will be using a planar approximation of 
+the poles. What this means is that we will map the 3D grid near the poles onto a 2D plane parallel to the poles, as if we put a giant flat plane in the exact center of the poles and draw lines
+from the grid directly upwards to the plane. For a visual representation, please consult the stream with timestamp 1:38:25 \cite{polarPlane}, which includes some explanation. In the streamm we
+use $r$ to indicate the radius of the planet (which we assume is a sphere), $\theta$ for the longitude and $\lambda$ for the latitude. So we have spherical coordinates, which we need to transform
+into $x$ and $y$ coordinates on the plane. We also need the distance between the center point (the point where the plane touches the planet which is the center of the pole) and the projected 
+point on the plane from the grid (the location on the plane where a line from the gird upwards to the plane hits it). This distance is denoted by $a$ (Simon chose this one, not me). We then get 
+the following equations as shown in \autoref{eq:polar distance}, \autoref{eq:polar x} and \autoref{polar y}. 
+
+\begin{subequations}
+    \begin{equation}
+        a = r \cos(\theta)
+        \label{eq:polar distance}
+    \end{equation}
+    \begin{equation}
+        x = a \sin(\lambda)
+        \label{eq:polar x}
+    \end{equation}
+    \begin{equation}
+        y = a \cos(\lambda)
+        \label{eq:polar y}
+    \end{equation}
+\end{subequations}
+
+But what if we know $x$ and $y$ and want to know $\theta$ and $\lambda$? Pythagoras' Theorem then comes into play \cite{pythagoras}. We know that (due to Pythagoras) \autoref{eq:pythagoras} must 
+always be true. Then if we substitue $a$ by $\sqrt{x^2 + y^2}$ in \autoref{eq:polar distance} we get \autoref{eq:polar theta1}. Then we transform that equation such that we only have $\theta$ on 
+one side and the rest on the other side (since we want to know $\theta$) and we get \autoref{eq:polar theta3}.
+\begin{equation}
+    x^2 + y^2 = a^2
+    \label{eq:pythagoras}
+\end{equation}
+
+\begin{subequations}
+    \begin{equation}
+        \sqrt{x^2 + y^2} = r\cos(\theta)
+        \label{eq:polar theta1}
+    \end{equation}
+    \begin{equation}
+        \frac{\sqrt{x^2 + y^2}}{r} = \cos(\theta)
+        \label{eq:polar theta2}
+    \end{equation}
+    \begin{equation}
+        \cos^{-1}(\frac{\sqrt{x^2 + y^2}}{r}) = \theta
+        \label{eq:polar theta3}
+    \end{equation}
+\end{subequations}
+
+For $\lambda$ we need another trigonometric function which is the tangent ($\tan$). The tangent is defined in \autoref{eq:tan}. If we then take a look at \autoref{eq:polar x} and 
+\autoref{eq:polar y}, we see that $\lambda$ is present in both equations. So we need to use both to get $\lambda$ \footnote{Yes you could only use one but since we both know $x$ and $y$ it is a
+bit easier to use both than to only use one as you need to know $\theta$ at that point as well which may or may not be the case.}. So let's combine \autoref{eq:polar x} and \autoref{eq:polar y}
+in \autoref{eq:polar lambda1}, transform it such that we end up with only $\lambda$ on one side and the rest on the other side and we end up with \autoref{eq:polar lambda3}.
+
+\begin{equation}
+    \tan(\alpha) = \frac{\sin(\alpha)}{\cos(\alpha)}
+    \label{eq:tan}
+\end{equation}
+
+\begin{subequations}
+    \begin{equation}
+        \frac{x}{y} = \frac{a\sin(\lambda)}{a\cos(\lambda)} = \frac{\sin(\lambda)}{\cos(\lambda)}
+        \label{eq:polar lambda1}
+    \end{equation}
+    \begin{equation}
+        \frac{x}{y} = \tan(\lambda)
+        \label{eq:polar lambda2}
+    \end{equation}
+    \begin{equation}
+        \lambda = \tan^{-1}(\frac{x}{y})
+        \label{eq:polar lambda3}
+    \end{equation}
+\end{subequations}
+
+With this math we can fix a lot of stuff in the model. With this we can resample (mapping from sphere to plane) the pressure, density, temperarature and advection to the plane and ensure that 
+there are no more overflows and funky business. The implementation (code) for this will be done in a follow up stream, so stay tuned!