# Moto Terra-'Luna' (in realtà la sua massa puà essere modificata),
# usando l'algoritmo che si trova su La Fisica di Feynman
# (Nota: i vettori sono pronti per varianti in 3D per prendere 
# in considerazione un terzo corpo (-> moti non coplanari)

source("dueCorpi_fun.R")

#-- Costanti --------------------------------------------------------
mT  = 5.97e24; mL  = 7.34e22    # kg
dTL = 384e6                     # m   (valore medio)
h   = 3600        # s   (1 ora)
vL0  = 1022       # m/s v nominale (https://it.wikipedia.org/wiki/Luna)

#-- Valori da variare per avere diversi moti
vL0 = 1.01*vL0    # la modifichiamo per avere diverse orbite 
#vL0 = 0.5*vL0    
# m = mL            # massa della Luna
# m = mT
m = mL/1000        # 'trascurabile' rispetto alla Terra 


dt = 1*h            #  discretizzazione
t.max = 12*30*24*h
iplot = 5           # ogni quanti step plottare i punti
CM = !TRUE          # per vedere il moto nel centro di massa del sistema

xlim=c(-1.1, 1.1)    # limiti del plot
ylim=c(-1.1, 1.1)

#--------------------------------------------------------------
# condizioni iniziali e velocità del centro di massa
# Terra
pT = c(0, 0, 0)       # vettori direttamente 3D in preparazione
vT = c(0, 0, 0)       # di eventuali n (ad es. n=3) corpi non coplanari
# Luna
pL = c(dTL, 0, 0)      
vL = c(0, vL0, 0)
# Centro di massa del sistema
vCM = (mT*vT + m*vL) / (mT + m)
#---------------------------------------------------------------

# iniziamo la simulazione
i = 1
t = 0
tv  = 0
dst = dist(pL, pT)

# plot vuoto, con un punto al centro
plot(0, 0, xlim=xlim,  ylim=ylim, 
     xlab='x / dTL', ylab='y / dTL', pch=4, cex=0.5,asp=1) 

# Forza SULLA Luna dovuta alla Terra
F = forza(m, pL, mT, pT)
cat(sprintf("F = (%e, %e)\n", F[1], F[2]))

# accelerazioni sui due corpi
aL =  F/m
aT = (-F)/mT    # terzo principio, reso esplicito mediante le parentesi     

cat("accelerazioni\n")
print(aL)
print(aT)

# calcoliamoci la velocità per t=-dt/2,
# al fine di usare lo stesso algoritmo di Feynman
vL = vL + aL * (-dt/2)
vT = vT + aT * (-dt/2)

cat("velocità\n")
print(vL)
print(vT)

cat("\nCominciamo: \n")

while(t < t.max) {
# velocità per t+dt/2 a partire a quelle per t-dt/2
  vL = vL + aL * dt
  vT = vT + aT * dt

  # posizioni per t+dt a partire da quelle per t
  pL = pL + vL * dt 
  pT = pT + vT * dt 

  if(i%%iplot == 0) { # plot ogni iplot passi per andare più veloci 
  # posizioni plottate in unità di distanze Terra-Luna  
  points( ( pL[1] + ifelse(CM, - vCM[1]*t, 0) ) / dTL,   
          ( pL[2] + ifelse(CM, - vCM[2]*t, 0) ) / dTL,
           col='orange', cex=0.5)
  points( ( pT[1] + ifelse(CM, - vCM[1]*t, 0) ) / dTL,
          ( pT[2] + ifelse(CM, - vCM[2]*t, 0) ) / dTL,
         col='blue', cex=0.5)
  }
  # forze e accelerazioni nelle nuove posizioni
  F = forza(m, pL, mT, pT)
  aL = F/m
  aT = (-F)/mT

  # incrementiamo il tempo
  t = t + dt
  i = i+1
  tv[i] = t
  dst[i] = dist(pL,pT)
}