cm: Completed exercise 1

3 files changed, 132 insertions, 59 deletions

diff --git a/Year 4/Celestial mechanics/Exercise1.R b/Year 4/Celestial mechanics/Exercise1.R
index 8e58c85..c76ce8a 100644
--- a/Year 4/Celestial mechanics/Exercise1.R
+++ b/Year 4/Celestial mechanics/Exercise1.R
@@ -1,14 +1,19 @@
+# Изменя величината спрямо подадено изменение и дата във векове от началото на епохата
 atDay <- function(epochConstant, rate, t) {
-  epochConstant + rate * t
+  epochConstant + t * rate
 }
 
+# Приема редица елементи и връща 3x3 матрица
 m3x3 <- function(...) {
   matrix(c(...), nrow = 3, ncol = 3, byrow = TRUE)
 }
+
+# Принтира аргументът, така че всички цифри да имат 8 знака прецизност
 p <- function(x) {
   formatC(x, format="f", digits = 8)
 }
 
+# Матрица на въртене
 Q <- function(theta, i, g) {
   m3x3(
     cos(theta), -sin(theta), 0,
@@ -25,6 +30,7 @@ Q <- function(theta, i, g) {
   )
 }
 
+# Векторът r
 coordinates <- function(theta, i, g, a, u, e) {
    Q(theta, i, g) %*% c(
     a * (cos(u) - e),
@@ -33,6 +39,7 @@ coordinates <- function(theta, i, g, a, u, e) {
   )
 }
 
+# Векторът v
 speed <- function(theta, i, g, a, u, e, n) {
   Q(theta, i, g) %*% c(
     -sin(u),
@@ -41,18 +48,22 @@ speed <- function(theta, i, g, a, u, e, n) {
   ) * 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) # защо конверитраме??
+  (dlambda - domega) / 100 / (2*pi)
 }
 
-
+# Ексцентрична аномалия
 u <- function(l, e, depth=8) {
   if (depth <= 0) {
     l
@@ -62,13 +73,14 @@ u <- function(l, e, depth=8) {
   }
 }
 
-planetData <- function(t,
-  a2000, da,
-  e2000, de,
-  i2000, di,
-  lambda2000, dlambda,
-  omega2000, domega,
-  theta2000, dtheta)
+# Дължина на вектор по Питагор
+length <- function(vec) {
+  sqrt(vec[1]^2 + vec[2]^2 + vec[3]^2)
+}
+
+# Координати, разстояние от слънцето, координатни скорости и скорост на планета по подадено
+# време и елиптични елементи
+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)
@@ -76,41 +88,48 @@ planetData <- function(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)
   
+  r <- coordinates(theta, i, g, a, u, e)
+  v <- speed(theta, i, g, a, u, e, n)
   
-  c(
-    coordinates(theta, i, g, a, u, e),
-    speed(theta, i, g, a, u, e, n)
-  )
+  c(r, length(r), v, length(v))
 }
 
-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)
-
-
-
+# Връща координати, разстояние от слънцето, координатни скорости и скорост на всяка планета
+# от слънчевата система, използвайки моета рожденна дата и J2000 данни
+solarSystemAtBirth <- function() {
+  birthDate <- 0.02320328
+  rbind(
+    planetData(birthDate,    # Меркурии
+               0.38709927,    0.00000037, 0.20563593, 0.00001906, 0.12225995, -0.00010380,
+               4.40259868, 2608.79030501, 1.35189358, 0.00280085, 0.84353100, -0.00218761),
+    planetData(birthDate,    # Венера
+               0.72333566,    0.00000390, 0.00677672, -0.00004107, 0.05924827, -0.00001377,
+               3.17613446, 1021.32854958, 2.29689636,  0.00004683, 1.33831572, -0.00484668),
+    planetData(birthDate,    # EM Bary
+               1.00000261,   0.00000562, 0.01671123, -0.00004392, -0.00000027, -0.00022596,
+               1.75343756, 628.30757790, 1.79660147,  0.00564219,  0.0,         0.0),
+    planetData(birthDate,    # Марс
+               1.52371034,   0.00001847,  0.09339410, 0.00007882, 0.03228321, -0.00014192,
+               -0.07947238, 334.06130168, -0.41789517, 0.00775643, 0.86497713, -0.00510637),
+    planetData(birthDate,    # Юпитер
+               5.20288700, -0.00011607, 0.04838624, -0.00013253, 0.02276602, -0.00003206,
+               0.60033114, 52.96631189, 0.25706047,  0.00370929, 1.75360053,  0.00357253),
+    planetData(birthDate,    # Сатурн
+               9.53667594, -0.00125060, 0.05386179, -0.00050991, 0.04338874,  0.00003379,
+               0.87186604, 21.33653879, 1.61615531, -0.00731244, 1.98378354, -0.00503838),
+    planetData(birthDate,    # Уран
+               19.18916464, -0.00196176, 0.04725744, -0.00004397, 0.01348507, -0.00004240,
+               5.46703627,  7.47842217, 2.98371499,  0.00712187, 1.29183904,  0.00074012),
+    planetData(birthDate,    # Нептун
+               30.06992276, 0.00026291, 0.00859048,  0.00005105, 0.03089309,  0.00000617,
+               -0.96202600, 3.81283674, 0.78478315, -0.00562720, 2.30006864, -0.00008878)
+  )
+}
 
+p(solarSystemAtBirth())
\ No newline at end of file
diff --git a/Year 4/Celestial mechanics/Project.pdf b/Year 4/Celestial mechanics/Project.pdf
index 2537d54..76cd281 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 10430d0..5800627 100644
--- a/Year 4/Celestial mechanics/Project.tex
+++ b/Year 4/Celestial mechanics/Project.tex
@@ -17,6 +17,7 @@
 \usepackage{amsfonts}           % Още математически символи
 \usepackage{cancel}             % Зачертаване
 \usepackage[dvipsnames]{xcolor} % Цветни кутийки + допълнителни цветове
+\usepackage{rotating}
 
 \usepackage{enumitem}
 \setlist[enumerate]{label=\arabic*.}
@@ -152,7 +153,7 @@ Neptune & 30.06992276   &   0.00859048   & 0.03089309  & -0.96202600   & 0.78478
 Дефинираме функция \texttt{atDay}, която пресмята новата стойност на величина спрямо поддадено $t$:
 \begin{lstlisting}[language=R]
 atDay <- function(epochConstant, rate, t) {
-  epochConstant + rate * t
+  epochConstant + t * rate
 }
 \end{lstlisting}
 
@@ -296,10 +297,12 @@ Eдин от начините да го дефинираме е чрез сре�
 	= \Delta \lambda - \Delta \tilde{\omega}
 \end{align*}
 
+Под тази формула, $n$ е в единици радиани за век, обаче на нас ни трябва безмерната ѝ стойност.
+Следователно, ще трябва да разделим $n$ на $100$ и $2 \pi$.
 Имплементираме формулата във функция \texttt{n}:
 \begin{lstlisting}[language=R]
 n <- function(dlambda, domega) {
-  (dlambda - domega) / 100 / (2*pi) # защо конверитраме??
+  (dlambda - domega) / 100 / (2*pi)
 }
 \end{lstlisting}
 
@@ -321,16 +324,22 @@ u <- function(l, e, depth=8) {
 }
 \end{lstlisting}
 
+Ще бъде нужно да пресметнем дължината на векторите $r$ и $v$.
+Дефинираме функция \texttt{length}, която приема 3-мерен вектор и връща дължината му, използвайки питагоровата формула:
+\begin{lstlisting}[language=R]
+length <- function(vec) { sqrt(vec[1]^2 + vec[2]^2 + vec[3]^2) }
+\end{lstlisting}
+
+\break
+
+\hfill
+
 Вече сме готови да имплементираме кода, който обединява всичко това.
-Дефинираме функция \texttt{planetData}, която приема елементите $a_{2000}$, $e_{2000}$, $i_{2000}$, $\lambda_{2000}$, $\tilde{\omega}_{2000}$, $\theta_{2000}$ и връща масив с две стойности: векторите $r$ и $v$:
+Дефинираме функция \texttt{planetData}, която приема датата $t$ във векове от началото на епохата и елементите $a_{2000}$, $e_{2000}$, $i_{2000}$, $\lambda_{2000}$, $\tilde{\omega}_{2000}$, $\theta_{2000}$ за една планета.
+Функцията пресмята координатите $x$, $y$, $z$; радиуса $|r|$; покомпонентните скорости $v_x$, $v_y$, $v_z$ и скоростта на планетата $v$:
+
 \begin{lstlisting}[language=R]
-planetData <- function(t,
-  a2000, da,
-  e2000, de,
-  i2000, di,
-  lambda2000, dlambda,
-  omega2000, domega,
-  theta2000, dtheta)
+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)
@@ -344,27 +353,72 @@ planetData <- function(t,
   n <- n(dlambda, domega)
   u <- u(l, e)
   
-  c(
-    coordinates(theta, i, g, a, u, e),
-    speed(theta, i, g, a, u, e, n)
-  )
+  r <- coordinates(theta, i, g, a, u, e)
+  v <- speed(theta, i, g, a, u, e, n)
+  
+  c(r, length(r), v, length(v))
 }
 \end{lstlisting}
 
+Финално дефинираме функцията \texttt{solarSystemAtBirth}, която извиква \texttt{planetData} за всяка планета с подадени константи:
+\begin{lstlisting}[language=R]
+solarSystemAtBirth <- function() {
+  birthDate <- 0.02320328
+  rbind(
+    planetData(birthDate,    # Меркурии
+               0.38709927,    0.00000037, 0.20563593, 0.00001906, 0.12225995, -0.00010380,
+               4.40259868, 2608.79030501, 1.35189358, 0.00280085, 0.84353100, -0.00218761),
+    planetData(birthDate,    # Венера
+               0.72333566,    0.00000390, 0.00677672, -0.00004107, 0.05924827, -0.00001377,
+               3.17613446, 1021.32854958, 2.29689636,  0.00004683, 1.33831572, -0.00484668),
+    planetData(birthDate,    # EM Bary
+               1.00000261,   0.00000562, 0.01671123, -0.00004392, -0.00000027, -0.00022596,
+               1.75343756, 628.30757790, 1.79660147,  0.00564219,  0.0,         0.0),
+    planetData(birthDate,    # Марс
+                1.52371034,   0.00001847,  0.09339410, 0.00007882, 0.03228321, -0.00014192,
+               -0.07947238, 334.06130168, -0.41789517, 0.00775643, 0.86497713, -0.00510637),
+    planetData(birthDate,    # Юпитер
+               5.20288700, -0.00011607, 0.04838624, -0.00013253, 0.02276602, -0.00003206,
+               0.60033114, 52.96631189, 0.25706047,  0.00370929, 1.75360053,  0.00357253),
+    planetData(birthDate,    # Сатурн
+               9.53667594, -0.00125060, 0.05386179, -0.00050991, 0.04338874,  0.00003379,
+               0.87186604, 21.33653879, 1.61615531, -0.00731244, 1.98378354, -0.00503838),
+    planetData(birthDate,    # Уран
+               19.18916464, -0.00196176, 0.04725744, -0.00004397, 0.01348507, -0.00004240,
+                5.46703627,  7.47842217, 2.98371499,  0.00712187, 1.29183904,  0.00074012),
+    planetData(birthDate,    # Нептун
+               30.06992276, 0.00026291, 0.00859048,  0.00005105, 0.03089309,  0.00000617,
+               -0.96202600, 3.81283674, 0.78478315, -0.00562720, 2.30006864, -0.00008878)
+  )
+}
+\end{lstlisting}
 
+Извикваме функцията:
+\begin{lstlisting}[language=R]
+p(solarSystemAtBirth())
+\end{lstlisting}
 
+И получаваме следния резултат:
 
-
-
-\begin{table}[!h]
-\begin{center}
+\begin{sidewaystable}
+\begin{center}\begin{Large}
 \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 \\
+Mercury & -0.25698453 & 0.21643004   & 0.04126659  & 0.33850545  & -1.38648452 & -1.18390004 & 0.03053503  & 1.82342837 \\
+Venus   & -0.11664062 & 0.70946840   & 0.01643289  & 0.71918043  & -1.16420814 & -0.19720358 & 0.06450187  & 1.18255247 \\
+EM Bary & -0.79848552 & -0.61311058  & 0.00000338  & 1.00671928  & 0.59280163  & -0.79699190 & 0.00000439  & 0.99328237 \\
+Mars    & 0.01018262  & 1.56865936   & 0.03261130  & 1.56903135  & -0.78255752 & 0.07431329  & 0.02078708  & 0.78635288 \\
+Jupiter & -1.81325044 & 4.88999123   & 0.02034819  & 5.21539121  & -0.41697241 & -0.13210905 & 0.00988286  & 0.43751167 \\
+Saturn  & 2.07612673  & 8.79599349   & -0.23579568 & 9.04076343  & -0.33293384 & 0.07398502  & 0.01195145  & 0.34126466 \\
+Uranus  & 16.49642718 & -11.30522661 & -0.25582952 & 20.00014268 & 0.12742638  & 0.17777408  & -0.00099186 & 0.21872835 \\
+Neptune & 18.92712374 & -23.39768921 & 0.04565204  & 30.09468321 & 0.14091098  & 0.11598482  & -0.00563534 & 0.18259282 \\
+
 \end{tabular}
-\end{center}
-\end{table}
+\end{Large}\end{center}
+\end{sidewaystable}
+
+\newpage
 
 \paragraph{2. Решение}