# Изменя величината спрямо подадено изменение и дата във векове от началото на епохата
atDay <- function(epochConstant, rate, t) {
  epochConstant + t * rate
}

# Принтира аргументът, така че всички цифри да имат 8 знака прецизност
p <- function(x) {
  formatC(x, format="f", digits = 21)
}

# Аргумент на перихелия
g <- function(omega, theta) {
  omega - theta
}

# Средна аномалия
l <- function(lambda, omega) {
  lambda - omega
}

gamma <- function(mu) { 1 + mu }
m <- function(mu) { mu / (1 + mu) }
L <- function(m, gamma, a) { m * sqrt(gamma) * sqrt(a) }
G <- function(L, e) { L * sqrt(1 - e^2) }
Theta <- function(G, i) { G * cos(i) }

# Елементи на Делоне по подадени константи
delaunay <- function(t, mu, 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)

  m <- m(mu)
  gamma <- gamma(mu)

  L <- L(m, gamma, a)
  G <- G(L, e)
  Theta <- Theta(G, i)
  l <- l(lambda, omega)
  g <- g(omega, theta)

  c(L, G, Theta, l, g, theta)
}

# Елементи на Делоне в моета дата на раждане
delaunayAtBirth <- function() {
  birthDate <- 0.023203285421
  rbind(
    delaunay(birthDate, 0.00000016601367952719,    # Меркурии
               0.38709927,    0.00000037, 0.20563593, 0.00001906, 0.12225995, -0.00010380,
               4.40259868, 2608.79030501, 1.35189358, 0.00280085, 0.84353100, -0.00218761),
    delaunay(birthDate, 0.00000244783833966454,    # Венера
               0.72333566,    0.00000390, 0.00677672, -0.00004107, 0.05924827, -0.00001377,
               3.17613446, 1021.32854958, 2.29689636,  0.00004683, 1.33831572, -0.00484668),
    delaunay(birthDate, 0.0000030404326333266,     # 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),
    delaunay(birthDate, 0.00000032271514450539,    # Марс
                1.52371034,   0.00001847,  0.09339410, 0.00007882, 0.03228321, -0.00014192,
               -0.07947238, 334.06130168, -0.41789517, 0.00775643, 0.86497713, -0.00510637),
    delaunay(birthDate, 0.00095479193842432214,    # Юпитер
               5.20288700, -0.00011607, 0.04838624, -0.00013253, 0.02276602, -0.00003206,
               0.60033114, 52.96631189, 0.25706047,  0.00370929, 1.75360053,  0.00357253),
    delaunay(birthDate, 0.0002858859806661029,     # Сатурн
               9.53667594, -0.00125060, 0.05386179, -0.00050991, 0.04338874,  0.00003379,
               0.87186604, 21.33653879, 1.61615531, -0.00731244, 1.98378354, -0.00503838),
    delaunay(birthDate, 0.00004366244043351564,    # Уран
               19.18916464, -0.00196176, 0.04725744, -0.00004397, 0.01348507, -0.00004240,
                5.46703627,  7.47842217, 2.98371499,  0.00712187, 1.29183904,  0.00074012),
    delaunay(birthDate, 0.00005151389020535497,    # Нептун
               30.06992276, 0.00026291, 0.00859048,  0.00005105, 0.03089309,  0.00000617,
               -0.96202600, 3.81283674, 0.78478315, -0.00562720, 2.30006864, -0.00008878)
  )
}

# Елементи на Поанкаре от първи род в моета дата на раждане
poincare1 <- function() {
  data <- delaunayAtBirth()
  result <- c()
  for (planet in 1:8) {
    L <- data[planet,1]; G <- data[planet,2]; Theta <- data[planet,3]
    l <- data[planet,4]; g <- data[planet,5]; theta <- data[planet,6]
    result <- rbind(
      result,
      c(L, L - G, G - Theta, l + g + theta, - g - theta, - theta)
    )
  }
  result
}

# Елементи на Поанкаре от втори род в моета дата на раждане
poincare2 <- function() {
  data <- delaunayAtBirth()
  result <- c()

  for (planet in 1:8) {
    L <- data[planet,1]; G <- data[planet,2]; Theta <- data[planet,3]
    l <- data[planet,4]; g <- data[planet,5]; theta <- data[planet,6]

    result <- rbind(
      result,
      c(L,              sqrt(2*(L-G)) * cos(g + theta),  sqrt(2*(G-Theta)) * cos(theta),
        l + g + theta, -sqrt(2*(L-G)) * sin(g + theta), -sqrt(2*(G-Theta)) * sin(theta))
    )
  }

  result
}

p(delaunayAtBirth())
p(poincare1())
p(poincare2())