cm: Mostly finished exercise 1

3 files changed, 379 insertions, 22 deletions

diff --git a/Year 4/Celestial mechanics/Exercise1.R b/Year 4/Celestial mechanics/Exercise1.R
new file mode 100644
index 0000000..8e58c85
--- /dev/null
+++ b/Year 4/Celestial mechanics/Exercise1.R
@@ -0,0 +1,116 @@
+atDay <- function(epochConstant, rate, t) {
+  epochConstant + rate * t
+}
+
+m3x3 <- function(...) {
+  matrix(c(...), nrow = 3, ncol = 3, byrow = TRUE)
+}
+p <- function(x) {
+  formatC(x, format="f", digits = 8)
+}
+
+Q <- function(theta, i, g) {
+  m3x3(
+    cos(theta), -sin(theta), 0,
+    sin(theta), cos(theta), 0,
+    0, 0, 1
+  ) %*% m3x3(
+    1, 0, 0,
+    0, cos(i), -sin(i),
+    0, sin(i), cos(i)
+  ) %*% m3x3(
+    cos(g), -sin(g), 0,
+    sin(g), cos(g), 0,
+    0, 0, 1
+  )
+}
+
+coordinates <- function(theta, i, g, a, u, e) {
+   Q(theta, i, g) %*% c(
+    a * (cos(u) - e),
+    a * sqrt(1 - e^2) * sin(u),
+    0
+  )
+}
+
+speed <- function(theta, i, g, a, u, e, n) {
+  Q(theta, i, g) %*% c(
+    -sin(u),
+    sqrt(1 - e^2) * cos(u),
+    0
+  ) * a * n / (1 - e * cos(u))
+}
+
+g <- function(omega, theta) {
+  omega - theta
+}
+l <- function(lambda, omega) {
+  lambda - omega
+}
+
+n <- function(dlambda, domega) {
+  (dlambda - domega) / 100 / (2*pi) # защо конверитраме??
+}
+
+
+u <- function(l, e, depth=8) {
+  if (depth <= 0) {
+    l
+  }
+  else {
+    l + e * sin(u(l, e, depth - 1))
+  }
+}
+
+planetData <- function(t,
+  a2000, da,
+  e2000, de,
+  i2000, di,
+  lambda2000, dlambda,
+  omega2000, domega,
+  theta2000, dtheta)
+{
+  a <- atDay(a2000, da, t)
+  e <- atDay(e2000, de, t)
+  i <- atDay(i2000, di, t)
+  lambda <- atDay(lambda2000, dlambda, t)
+  omega <- atDay(omega2000, domega, t)
+  theta <- atDay(theta2000, dtheta, t)
+  print(a)
+  print(e)
+  print(i)
+  print(lambda)
+  print(omega)
+  print(theta)
+  
+  g <- g(omega, theta)
+  l <- l(lambda, omega)
+  n <- n(dlambda, domega)
+  u <- u(l, e)
+  
+  
+  c(
+    coordinates(theta, i, g, a, u, e),
+    speed(theta, i, g, a, u, e, n)
+  )
+}
+
+p(planetData(0.02320328,
+           1.00000261, 0.00000562,
+           0.01671123, -0.00004392,
+           -0.00000027, -0.00022596,
+           1.75343756, 628.30757790,
+           1.79660147, 0.00564219,
+           0,0))
+
+planetData(0.02320328,
+             0.38709927, 0.00000037,
+             0.20563593, 0.00001906,
+             0.12225995, -0.00010380,
+             4.40259868, 2608.79030501,
+             1.35189358, 0.00280085,
+             0.84353100, -0.00218761)
+
+
+
+
diff --git a/Year 4/Celestial mechanics/Project.pdf b/Year 4/Celestial mechanics/Project.pdf
index 8c0c8ce..2537d54 100644
--- a/Year 4/Celestial mechanics/Project.pdf
+++ b/Year 4/Celestial mechanics/Project.pdf
diff --git a/Year 4/Celestial mechanics/Project.tex b/Year 4/Celestial mechanics/Project.tex
index b8f671a..10430d0 100644
--- a/Year 4/Celestial mechanics/Project.tex
+++ b/Year 4/Celestial mechanics/Project.tex
@@ -21,6 +21,25 @@
 \usepackage{enumitem}
 \setlist[enumerate]{label=\arabic*.}
 
+\usepackage{listings}
+\newfontfamily{\lstsansserif}[Scale=1]{DejaVuSansMono}
+\lstdefinestyle{mystyle}{
+    basicstyle=\lstsansserif\footnotesize,
+    breakatwhitespace=false,
+    breaklines=true,
+    captionpos=b,
+    keepspaces=true,
+    numbers=left,
+    numbersep=15pt,
+    showspaces=false,
+    showstringspaces=false,
+    showtabs=false,
+    tabsize=2,
+    frame=single,
+    literate = {-}{-}1{"}{\textquotedbl}1{*}{\normalfont{*}}1,
+}
+\lstset{style=mystyle}
+
 \newcommand{\rad}{\mathrm{rad}}
 
 \usepackage[pdfauthor={Камен Димитров Младенов, },
@@ -74,57 +93,279 @@ Uranus  & 19.18916464   &   0.04725744   &   0.77263783  &    313.23810451  &  1
 \hline
 Neptune & 30.06992276   &   0.00859048   &   1.77004347  &    -55.12002969  &   44.96476227  &  131.78422574 \\
         &  0.00026291   &   0.00005105   &   0.00035372  &    218.45945325  &   -0.32241464  &   -0.00508664 \\
-
 \end{tabular}
 \end{center}
 \end{table}
 
 Първо, всички колони със стойности спрямо градуси ще бъдат преобразувани спрямо радиани, използвайки формулата $\rad = \deg \cdot \frac{\pi}{180}$ (с точност до 8 знака след десетичния знак).
 
-Второ, ще сменим наименованието на някои константи: $\mathrm{I} \rightarrow \mathrm{i}$, $\mathrm{L} \rightarrow \mathrm{\lambda}$, $\mathrm{long.peri.} \rightarrow \tilde{\omega}$ и \\
-$\mathrm{long.node.} \rightarrow \mathrm{\theta}$.
-Получаваме следната таблица:
+Второ, ще сменим наименованиетa на елементите: $\mathrm{a} \rightarrow a_{2000}$, $\mathrm{e} \rightarrow e_{2000}$, $\mathrm{I} \rightarrow i_{2000}$, $\mathrm{L} \rightarrow \lambda_{2000}$, $\mathrm{long.peri.} \rightarrow \tilde{\omega}_{2000}$ и $\mathrm{long.node.} \rightarrow \theta_{2000}$.
+Ще бележим измененията с $\Delta$.
+Получаваме:
 
 \begin{table}[!h]
 \begin{center}
 \begin{tabular}{c||c|c|c|c|c|c} 
-        &       a       &       e        &       i       &         $\lambda$        &    $\tilde{\omega}$  &    $\theta$ \\
+        &       $a_{2000}, \Delta a$       &       $e_{2000}, \Delta e$        &       $i_{2000}, \Delta i$       &         $\lambda_{2000}, \Delta \lambda$        &    $\tilde{\omega}_{2000}, \Delta \tilde{\omega}$  &    $\theta_{2000}, \Delta \theta$ \\
         &   au, au/Cy   &  rad, rad/Cy   &  rad, rad/Cy  &    rad, rad/Cy   &   rad, rad/Cy  &   rad, rad/Cy \\
 \hline\hline
-Mercury &  0,38709927   &   0,20563593   & 0.12225995  & 4.40259868    & 1.35189358  & 0.84353100 \\
-        &  0,00000037   &   0,00001906   & -0.00010380 & 2608.79030501 & 0.00280085  & -0.00218761 \\
+Mercury &  0.38709927   &   0.20563593   & 0.12225995  & 4.40259868    & 1.35189358  & 0.84353100 \\
+        &  0.00000037   &   0.00001906   & -0.00010380 & 2608.79030501 & 0.00280085  & -0.00218761 \\
 \hline
-Venus   &  0,72333566   &   0,00677672   & 0.05924827  & 3.17613446    & 2.29689636  & 1.33831572 \\
-        &  0,00000390   &  -0,00004107   & -0.00001377 & 1021.32854958 & 0.00004683  & -0.00484668 \\
+Venus   &  0.72333566   &   0.00677672   & 0.05924827  & 3.17613446    & 2.29689636  & 1.33831572 \\
+        &  0.00000390   &  -0.00004107   & -0.00001377 & 1021.32854958 & 0.00004683  & -0.00484668 \\
 \hline
-EM Bary &  1,00000261   &   0,01671123   & -0.00000027 & 1.75343756    & 1.79660147  & 0.0        \\
-        &  0,00000562   &  -0,00004392   & -0.00022596 & 628.30757790  & 0.00564219  & 0.0        \\
+EM Bary &  1.00000261   &   0.01671123   & -0.00000027 & 1.75343756    & 1.79660147  & 0.0        \\
+        &  0.00000562   &  -0.00004392   & -0.00022596 & 628.30757790  & 0.00564219  & 0.0        \\
 \hline
-Mars    &  1,52371034   &   0,09339410   & 0.03228321  & -0.07947238   & -0.41789517 & 0.86497713 \\
-        &  0,00001847   &   0,00007882   & -0.00014192 & 334.06130168  & 0.00775643  & -0.00510637 \\
+Mars    &  1.52371034   &   0.09339410   & 0.03228321  & -0.07947238   & -0.41789517 & 0.86497713 \\
+        &  0.00001847   &   0.00007882   & -0.00014192 & 334.06130168  & 0.00775643  & -0.00510637 \\
 \hline
-Jupiter &  5,20288700   &   0,04838624   & 0.02276602  & 0.60033114    & 0.25706047  & 1.75360053 \\
-        & -0,00011607   &  -0,00013253   & -0.00003206 & 52.96631189   & 0.00370929  & 0.00357253 \\
+Jupiter &  5.20288700   &   0.04838624   & 0.02276602  & 0.60033114    & 0.25706047  & 1.75360053 \\
+        & -0.00011607   &  -0.00013253   & -0.00003206 & 52.96631189   & 0.00370929  & 0.00357253 \\
 \hline
-Saturn  &  9,53667594   &   0,05386179   & 0.04338874  & 0.87186604    & 1.61615531  & 1.98378354 \\
-        & -0,00125060   &  -0,00050991   & 0.00003379  & 21.33653879   & -0.00731244 & -0.00503838 \\
+Saturn  &  9.53667594   &   0.05386179   & 0.04338874  & 0.87186604    & 1.61615531  & 1.98378354 \\
+        & -0.00125060   &  -0.00050991   & 0.00003379  & 21.33653879   & -0.00731244 & -0.00503838 \\
 \hline
-Uranus  & 19,18916464   &   0,04725744   & 0.01348507  & 5.46703627    & 2.98371499  & 1.29183904 \\
-        & -0,00196176   &  -0,00004397   & -0.00004240 & 7.47842217    & 0.00712187  & 0.00074012 \\
+Uranus  & 19.18916464   &   0.04725744   & 0.01348507  & 5.46703627    & 2.98371499  & 1.29183904 \\
+        & -0.00196176   &  -0.00004397   & -0.00004240 & 7.47842217    & 0.00712187  & 0.00074012 \\
 \hline
-Neptune & 30,06992276   &   0,00859048   & 0.03089309  & -0.96202600   & 0.78478315  & 2.30006864 \\
-        &  0,00026291   &   0,00005105   & 0.00000617  & 3.81283674    & -0.00562720 & -0.00008878 \\
+Neptune & 30.06992276   &   0.00859048   & 0.03089309  & -0.96202600   & 0.78478315  & 2.30006864 \\
+        &  0.00026291   &   0.00005105   & 0.00000617  & 3.81283674    & -0.00562720 & -0.00008878 \\
 \end{tabular}
 \end{center}
 \end{table}
 
-Също като подготовка, ще изчислим рожденната ми дата във векове.
+Без долен индекс и $\Delta$ ще обозначаваме стойността на съответната величина след време $t$ от началото на епохата.
+Тоест, обозначаваме $a = a_{2000} + t \cdot \Delta a$, $e = e_{2000} + t \cdot \Delta e$, ...
+
+Финално, ще изчислим рожденната ми дата във векове.
 Аз съм роден на 27 Април 2002 година.
 Това е 117-ят ден в годината, което е $\frac{117}{365.25} = 0.32032854$ Юлиански години.
 Тоест, роден съм $2.32032854$ Юлиански години след началото на епохата (1 Януари 2000, 00:00 часа), което са $\frac{2.32032854}{100} = 0.02320328$ века.
 
 \paragraph{1. Решение}
 
+Ще напишем програма на език \texttt{R}, която намира координатите и скоростите на планетите в деня на моето раждане.
+Единственият вход на програмата ще бъде горепосочената таблица и рожденната ми дата (в Юлиански векове от началото на епохата).
+
+Дефинираме функция \texttt{atDay}, която пресмята новата стойност на величина спрямо поддадено $t$:
+\begin{lstlisting}[language=R]
+atDay <- function(epochConstant, rate, t) {
+  epochConstant + rate * t
+}
+\end{lstlisting}
+
+От теорема знаем, че векторът координат $r$ и векторът скорост $v$ се пресмятат чрез:
+\begin{align*}
+	r = \begin{pmatrix}
+		x \\ y \\ z \\
+	\end{pmatrix}
+	= Q
+	\begin{pmatrix}
+		a(\cos u - e) \\
+		a \sqrt{1 - e^2} \sin u \\
+		0
+	\end{pmatrix}
+	\qquad
+	v = \begin{pmatrix}
+		v_x \\ v_y \\ v_z \\
+	\end{pmatrix}
+	= Q
+	\begin{pmatrix}
+		-\sin u \\
+		\sqrt{1 - e^2} \cos u \\
+		0
+	\end{pmatrix}
+	\frac{an}{1 - e \cos u}
+\end{align*}
+
+Където $Q$ е матрицата на въртения около нулата:
+\begin{align*}
+	Q
+	= \begin{pmatrix}
+		\cos\theta & -\sin\theta & 0 \\
+		\sin\theta & \cos\theta & 0 \\
+		0 & 0 & 1 \\
+	\end{pmatrix}
+	\begin{pmatrix}
+		1 & 0 & 0 \\
+		0 & \cos i & -\sin i \\
+		0 & \sin i & \cos i \\
+	\end{pmatrix}
+	\begin{pmatrix}
+		\cos g & -\sin g & 0 \\
+		\sin g & \cos g & 0 \\
+		0 & 0 & 1 \\
+	\end{pmatrix}
+\end{align*}
+
+Дефинираме помощни функции \texttt{m3x3} и \texttt{p}. Първата приема редица елементи и ги преобразува в $3 \times 3$ матрица, втората изкарва подаденото на екрана, като всички числа са с 8 цифри прецизност:
+
+\begin{lstlisting}[language=R]
+m3x3 <- function(...) {
+  matrix(c(...), nrow = 3, ncol = 3, byrow = TRUE)
+}
+p <- function(x) {
+  formatC(x, format="f", digits = 8)
+}
+\end{lstlisting}
+
+Дефинираме функция \texttt{Q}, която пресмята матрица на въртене по подадени стойности $\theta$, $i$, $g$ (операторът \texttt{\%*\%} умножава матрици и вектори):
+
+\begin{lstlisting}[language=R]
+Q <- function(theta, i, g) {
+  m3x3(
+    cos(theta), -sin(theta), 0,
+    sin(theta),  cos(theta), 0,
+             0,           0, 1
+  ) %*% m3x3(
+    1,      0,       0,
+    0, cos(i), -sin(i),
+    0, sin(i),  cos(i)
+  ) %*% m3x3(
+    cos(g), -sin(g), 0,
+    sin(g),  cos(g), 0,
+         0,       0, 1
+  )
+}
+\end{lstlisting}
+
+Дефинираме функция \texttt{coordinates}, която изчислява стойността на $r$ вектора:
+
+\begin{lstlisting}[language=R]
+coordinates <- function(theta, i, g, a, u, e) {
+   Q(theta, i, g) %*% c(
+    a * (cos(u) - e),
+    a * sqrt(1 - e^2) * sin(u),
+    0
+  )
+}
+\end{lstlisting}
+
+Дефинираме функция \texttt{speed}, която изчислява стойността на $v$ вектора:
+
+\begin{lstlisting}[language=R]
+speed <- function(theta, i, g, a, u, e, n) {
+   Q(theta, i, g) %*% c(
+    -sin(u),
+    sqrt(1 - e^2) * cos(u),
+    0
+  ) * a * n / (1 - e * cos(u))
+}
+\end{lstlisting}
+
+По дефиниция знаем, че $\tilde{\omega} = g + \theta$ и $\lambda = l + g + \theta = l + \tilde{\omega}$.
+То, можем да намерим стойностите на $g$ и $l$ чрез формулите:
+\begin{align*}
+	g = \tilde{\omega} - \theta
+	\qquad
+	l = \lambda - \tilde{\omega}
+\end{align*}
+
+Имплементираме ги във функции \texttt{g} и \texttt{l}:
+\begin{lstlisting}[language=R]
+g <- function(omega, theta) { omega - theta }
+l <- function(lambda, omega) { lambda - omega }
+\end{lstlisting}
+
+Величанта $n$ е "средното движение".
+Eдин от начините да го дефинираме е чрез средната аномалия:
+\begin{align*}
+	l = l_0 + n (t - t_0)
+	\qquad
+	n = \frac{l - l_0}{t - t_0}
+\end{align*}
+
+В нашата задача, $l_0 = l_{2000} = \lambda_{2000} - \tilde{\omega}_{2000}$.
+Нека да зачетем $t_0$ за началото на епохата и $t$ за време от това начало (във векове), тоест свеждаме $t - t_0$ до $t$.
+Тогава намираме, че $n$ е изменението на средната аномалия $l$:
+\begin{align*}
+	n = \frac{l - l_{2000}}{t}
+	= \frac{(l_{2000} + t \cdot \Delta l) - l_{2000}}{t}
+	= \frac{t \cdot \Delta l}{t}
+	= \Delta l
+\end{align*}
+
+Ще покажем, че изменението на средната аномалия е измението на $\lambda$ минус изменението на $\tilde{\omega}$:
+\begin{align*}
+	\Delta l & = n = \frac{l - l_{2000}}{t}
+	= \frac{(\lambda - \tilde{\omega}) - (\lambda_{2000} - \tilde{\omega}_{2000})}{t}
+	= \frac{(\lambda_{2000} + t \cdot \Delta \lambda) - (\tilde{\omega}_{2000} + t \cdot \Delta \tilde{\omega}) - \lambda_{2000} + \tilde{\omega}_{2000}}{t} \\
+	& = \frac{t \cdot \Delta \lambda - t \cdot \Delta \tilde{\omega}}{t}
+	= \Delta \lambda - \Delta \tilde{\omega}
+\end{align*}
+
+Имплементираме формулата във функция \texttt{n}:
+\begin{lstlisting}[language=R]
+n <- function(dlambda, domega) {
+  (dlambda - domega) / 100 / (2*pi) # защо конверитраме??
+}
+\end{lstlisting}
+
+Ще дефинираме стойността $u$ използвайки вторият метод за решаването на уравнението на Кеплер:
+\begin{align*}
+	u = l + e \cdot \sin(l + e \cdot \sin (l + e \cdot \sin( \cdots )))
+\end{align*}
+
+Нуждаем се от 8 знака прецизност, затова ще използваме следното приближение:
+\begin{align*}
+	u = l + e \cdot \sin(l + e \cdot \sin (l + e \cdot \sin(l + e \cdot \sin(l + e \cdot \sin(l + e \cdot \sin(l + e \cdot \sin(l + e \cdot \sin(l))))))))
+\end{align*}
+
+Дефинираме функция \texttt{u}, която пресмята стойността на $u$ чрез рекурсия:
+\begin{lstlisting}[language=R]
+u <- function(l, e, depth=8) {
+  if (depth <= 0) { l }
+  else            { l + e * sin(u(l, e, depth - 1)) }
+}
+\end{lstlisting}
+
+Вече сме готови да имплементираме кода, който обединява всичко това.
+Дефинираме функция \texttt{planetData}, която приема елементите $a_{2000}$, $e_{2000}$, $i_{2000}$, $\lambda_{2000}$, $\tilde{\omega}_{2000}$, $\theta_{2000}$ и връща масив с две стойности: векторите $r$ и $v$:
+\begin{lstlisting}[language=R]
+planetData <- function(t,
+  a2000, da,
+  e2000, de,
+  i2000, di,
+  lambda2000, dlambda,
+  omega2000, domega,
+  theta2000, dtheta)
+{
+  a <- atDay(a2000, da, t)
+  e <- atDay(e2000, de, t)
+  i <- atDay(i2000, di, t)
+  lambda <- atDay(lambda2000, dlambda, t)
+  omega <- atDay(omega2000, domega, t)
+  theta <- atDay(theta2000, dtheta, t)
+  
+  g <- g(omega, theta)
+  l <- l(lambda, omega)
+  n <- n(dlambda, domega)
+  u <- u(l, e)
+  
+  c(
+    coordinates(theta, i, g, a, u, e),
+    speed(theta, i, g, a, u, e, n)
+  )
+}
+\end{lstlisting}
+
+
+
+
+
+
+\begin{table}[!h]
+\begin{center}
+\begin{tabular}{c||c|c|c|c|c|c|c|c} 
+& $x$ & $y$ & $z$ & $|r|$ & $v_x$ & $v_y$ & $v_z$ & $|v|$ \\
+\hline\hline
+Mercury \\
+\end{tabular}
+\end{center}
+\end{table}
+
 \paragraph{2. Решение}
 
 \end{document}
\ No newline at end of file