# Trasformazione del sistema Terra-Luna in un 'tre corpi'

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 (https://it.wikipedia.org/wiki/Luna)

#-- Valori da variare per avere diversi moti
# corpo 3 è quello che ha velocità diversa da zero
# come la Luna nel caso di due_corpi.R

m1 = mT

m2 = mT/100

m3 = mT/100

dt = 1*h            #  discretizzazione
t.max = 12*30*24*h
iplot = 10           # ogni quanti passi fi plottano i punti
CM = TRUE           # per vedere il moto nel centro di massa del sistema

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

#--------------------------------------------------------------
# condizioni iniziali e velocità del centro di massa
# Terra e altro corpo
p1 = c(0, 0, 0)
v1 = c(0, 0, 0)
p2 = c(-dTL, 0, 0)
v2 = c(0, -vL0/2, 0)

# 'Luna'
p3 = c(dTL, 0, 0)      
v3 = c(0, vL0, 0)

# Centro di massa del sistema
vCM = (m1*v1 + m2*v2 + m3*v3) / (m1 + m2 + m3)
#---------------------------------------------------------------

# iniziamo la simulazione
i = 1
t = 0

# 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 sui ciascun corpo dovuto agli altri due
F1 = forza(m1, p1, m2, p2) + forza(m1, p1, m3, p3)
F2 = forza(m2, p2, m1, p1) + forza(m2, p2, m3, p3)
F3 = forza(m3, p3, m1, p1) + forza(m3, p3, m2, p2)

# accelerazioni sui tre corpi
a1 =  F1/m1
a2 =  F2/m2
a3 =  F3/m3

# calcoliamoci la velocità per t=-dt/2,
# al fine di usare lo stesso algoritmo di Feynman
v1 = v1 + a1 * (-dt/2)
v2 = v2 + a2 * (-dt/2)
v3 = v3 + a3 * (-dt/2)


cat("\nCominciamo: \n")

while(t < t.max) {
  # velocità per t+dt/2 a partire a quelle per t-dt/2
  v1 = v1 + a1 * dt
  v2 = v2 + a2 * dt    
  v3 = v3 + a3 * dt
    
  # posizioni per t+dt a partire da quelle per t
  p1 = p1 + v1 * dt 
  p2 = p2 + v2 * dt
  p3 = p3 + v3 * dt

  if(i%%iplot == 0) { # plot ogni iplot passi per andare più veloci 
     # posizioni plottate in unità di distanze Terra-Luna  
     points( ( p1[1] + ifelse(CM, - vCM[1]*t, 0) ) / dTL,
             ( p1[2] + ifelse(CM, - vCM[2]*t, 0) ) / dTL,
            col='blue', cex=0.5)
     points( ( p2[1] + ifelse(CM, - vCM[1]*t, 0) ) / dTL,
             ( p2[2] + ifelse(CM, - vCM[2]*t, 0) ) / dTL,
            col='red', cex=0.5)
     points( ( p3[1] + ifelse(CM, - vCM[1]*t, 0) ) / dTL,
             ( p3[2] + ifelse(CM, - vCM[2]*t, 0) ) / dTL,
            col='orange', cex=0.5)
  }
  # forze e accelerazioni nelle nuove posizioni
  F1 = forza(m1, p1, m2, p2) + forza(m1, p1, m3, p3)
  F2 = forza(m2, p2, m1, p1) + forza(m2, p2, m3, p3)
  F3 = forza(m3, p3, m1, p1) + forza(m3, p3, m2, p2)
  a1 =  F1/m1
  a2 =  F2/m2
  a3 =  F3/m3

  # incrementiamo il tempo
  t = t + dt
  i = i+1
}