#---------------------------------------------------------
# Precessione (*) di momento della quantità di moto
# e di momento della forza (ruota-giroscopio)   
#
# GdA, maggio 2024
#--------------------------------------------------------

# funzione che disegna un cerchio 
cerchio <- function(C=c(0,0), r=1, col='blue') {
  th <- seq(0, 2*pi, by=0.02)
  points(C[1] + r*cos(th), C[2] + r*sin(th), ty='l', col)
}

# plot vuoto
xyM <- 10
plot(NULL,  xaxt = 'n', yaxt = 'n', bty = 'n', pch = '', asp=1,
     ylab = '', xlab = '', xlim=c(-xyM,xyM), ylim=c(-xyM,xyM))

# disegniamo il cerchio
R <- 3
cerchio(r=R, col='gray')

# decidiamo il periodo (*) e tracciamo le frecce a ogmi frazione di periodo 
T0  <- 5   # s
# T0  <- 5/2   # un periodo minore corrisponde, a parità di L, a M maggiore!
ndt <- 16
dt  <- T0/ndt
om0 <- 2*pi/T0    # (*)
sM <- 1; sL <- 3  # fattori di scala per plottare i vettori M  e L
lwd = 1.2
for (i in 0:(2*ndt)) {    # facciamo due giri     
    if (i==ndt) lwd=2*lwd # al secondo giro intensifichiamo le frecce
    t <- i*dt
    r <- R * c( cos(om0*t), sin(om0*t) )    # punti sulla circonferenza
    L <- sL * c( cos(om0*t), sin(om0*t) )   # vettore L
    points(r[1],r[2], pch=19, cex=0.5)
    arrows(r[1], r[2], r[1]+L[1], r[2]+L[2], len=0.1, col='blue' , lwd=lwd)
    
    M <- om0*R* c(-sin(om0*t), cos(om0*t))  # vettore M 
    arrows(r[1], r[2], r[1]+sM*M[1], r[2]+sM*M[2], len=0.1, col='red' , lwd=lwd)
    Sys.sleep(dt)  # pausa di 'dt' 
}

#------------------------------------------------------------------------------
# Nota: T0 e quindi om0 si riferiscono alla precessione dei vettori L e M 
# e NON alla rotazione della ruota-giroscopio intorno al proprio asse!
# ------------------------------------------------------------------------------