# Изменя величината спрямо подадено изменение и дата във векове от началото на епохата
atDay <- function(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 = 10)
}

# Матрица на въртене
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
  )
}

# Векторът 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
  )
}

# Векторът v
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))
  }
}

# Дължина на вектор по Питагор
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)
  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)
  
  r <- coordinates(theta, i, g, a, u, e)
  v <- speed(theta, i, g, a, u, e, n)
  
  c(r, length(r), v, length(v))
}

# Връща координати, разстояние от слънцето, координатни скорости и скорост на всяка планета
# от слънчевата система, използвайки моета рожденна дата и J2000 данни
solarSystemAtBirth <- function() {
  birthDate <- 0.023203285421
  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())