# angolo.fra.vettori <- function(v1, v2) {
#  # angolo fra v1 e v2 (nel piano)
#  # calcolato dal prodotto scalare diviso i moduli dei vettori
#  acos(   as.numeric( t(v1) %*% v2 )  /
#        sqrt( as.numeric( t(v1) %*% v1 )  * as.numeric( t(v2) %*% v2 ) )
#      ) * 180/pi 
# }

# VERSIONE AGGIORNATA, meno criptica

# prodotto scalare:
sc.prod <- function(v1, v2)  as.numeric( t(v1) %*% v2 )

# angolo fra vettori
angolo.fra.vettori <- function(v1, v2) {
  costh <- sc.prod(v1, v2) / (sqrt(sc.prod(v1, v1)) * sqrt(sc.prod(v2, v2)))
  return( acos(costh) * 180/pi )
}

#--------------------------------

R.T = 1 # U.A    (unità astronomica)
T.T = 1 # anno

R.G = 5.2 #
T.G = T.T * R.G^(3/2)

omega.T = 2*pi/T.T
omega.G = 2*pi/T.G

cat(sprintf("Distanza Giove-Sole: %.2fU.A; periodo: %.2f y\n", R.G, T.G))

cat(sprintf("Vel. angolari: Terra: %.3f rad/y;  Giove: %.2f rad/y\n",
            omega.T, omega.G))

t <- seq(0, 1, by=1/104)

xT = R.T * cos(omega.T*t)
yT = R.T * sin(omega.T*t)
plot(xT, yT, asp=1, col='blue')
readline(prompt = "Premi Invio per continuare")

t <- seq(0, 12, by=1/104)
xG = R.G * cos(omega.G*t)
yG = R.G * sin(omega.G*t)

xT = R.T * cos(omega.T*t)
yT = R.T * sin(omega.T*t)

plot(xG, yG, asp=1, col='orange')
points(xT, yT, cex=0.5, col='blue')

readline(prompt = "Premi Invio per continuare")
plot(xG-xT, yG-yT, asp=1, col='orange')

readline(prompt = "Premi Invio per continuare")
plot(xG[1]-xT[1], yG[1]-yT[1], asp=1, col='orange',
    xlab='xG - xT', ylab='yG - yT',
    xlim=range(xG-xT), ylim=range(yG-yT))  
for(i in 2:length(t)) {
  points(xG[i]-xT[i], yG[i]-yT[i], col='orange')
  Sys.sleep(0.002)
}

# Complemento: plot solo se Giove è 'ben visibile'
readline(prompt = "Premi Invio per continuare")
plot(xG[1]-xT[1], yG[1]-yT[1], asp=1, col='orange',
    xlab='xG - xT', ylab='yG - yT',
    xlim=range(xG-xT), ylim=range(yG-yT))  
for(i in 2:length(t)) {
  if (angolo.fra.vettori(c(xT[i],yT[i]), c(xG[i],yG[i])) < 45) { 
     colore = 'orange'    
  } else {
     colore = 'gray97'   # se non visibile
  }
  points(xG[i]-xT[i], yG[i]-yT[i], col=colore)
  Sys.sleep(0.002)
}