#---------------------------------------------------
# Oscillazione della molla: valutazione numerica
#
# Versione con 'animazione'
#
# GdA, marzo 2025
#---------------------------------------------------

k  = 50      # (N/m)   costante della molla
m  = 0.5     # kg      massa sospesa
xM = 0.05    # m       allungamento massimo (a t=0)

dt   = 0.01  # s       intervallo di discretizzazione
tMax = 1.5   # s       tempo massimo 

t = 0                                                      # 1)
x = xM
v = 0
a = - k * x / m

i = 2  # indice delle iterazioni

while( t[i-1] < tMax ) {                                   
   t[i] <- t[i-1] + dt                                     
   v.m0 <- v[i-1] + a[i-1]*dt/2                            # 2)
   x.m0 <- x[i-1] + v.m0*dt/2                              # 3)
   a.m  <- ( - k * x.m0 ) / m                              # 4) 
   v[i] <- v[i-1] + a.m*dt                                 # 5)
   x[i] <- x[i-1] + (v[i-1] + v[i])/2 * dt                 # 6)
   a[i] <- ( - k * x[i] ) / m                              # 7)
   i <- i + 1
}

plot(NA, NA, pch=19, col='blue', cex=0.5, xaxs = "i", yaxs = "i",
     xlab='t (s)', ylab='x/xMax;  v/vMax;  a/aMax',
     xlim=c(0, tMax), ylim=c(-1,1.2) )
text(0.15*max(t), 1.15, sprintf("xMax = %.2f m", max(x)), col='blue', cex=1.7)
text(0.5*max(t), 1.15, sprintf("vMax = %.2f m/s", max(v)), col='cyan', cex=1.7)
text(0.85*max(t), 1.15, sprintf("aMax = %.2f m/s^2", max(a)), col='red', cex=1.7)
abline(h=0, col='gray')
text(0.5*max(t), 1.06, sprintf("Soluzione numerica"), col='black', cex=2.5)
# T=2*pi/sqrt(k/m)
# abline(v=T)
# abline(v=2*T)
xMax <- max(x)
vMax <- max(v)
aMax <- max(a)
for(i in 1:length(t)) {
  points(t[i], x[i]/xMax,  pch=19, col='blue', cex=0.5)  
  points(t[i], v[i]/vMax,  pch=19, col='cyan', cex=0.5)
  points(t[i], a[i]/aMax,  pch=19, col='red', cex=0.5)
  Sys.sleep(dt)
}
#Sys.sleep(1)
#points(t, x/xMax,  ty='l', lty=2, col='blue') 
#Sys.sleep(1)
#points(t, v/vMax,  ty='l', lty=2, col='cyan') 
#Sys.sleep(1)
#points(t, a/aMax,  ty='l', lty=2, col='red')

ti <- seq(0, tMax, len=201)
Sys.sleep(1.5)
text(0.5*max(t), 1.06, sprintf("Soluzione numerica"), col='white', cex=2.5)
text(0.5*max(t), 1.06, sprintf("Soluzione esatta"), col='black', cex=2.5)
Sys.sleep(0.5)
omega = sqrt(k/m)
points(ti, cos(omega*ti),  ty='l', lty=1, lwd=2, col='blue') 
Sys.sleep(1)
points(ti, -sin(omega*ti),  ty='l', lty=1, lwd=2, col='cyan') 
Sys.sleep(1)
points(ti, -cos(omega*ti),  ty='l', lty=1, lwd=2, col='red')
#-- Note ------------------

# 1) Si ricorda che in R le variabili sono per default dei vettori,
#   ovvero "t = 0" crea il vettore t e assegna 0 al primo elemento.
#   Quindi prima di entrare nel loop while t e t[1] sono la stessa cosa
#   Idem per x, v e a
#
# 2) Si calcola la velocità a metà del nuovo intervallino dt 
#    usando come accelerazione quella calcolata alla posizione
#    del tempo precedente
#     Nota: v[i-1] + a[i-1]*dt/2 = (v[i-1] + v[i-1] + a[i-1]*dt)/2 
#
# 3) Posizione a metà del nuovo intervallino, basata sulla velocità
#    'intermedia' v.m0
#     (Vale nota analoga al punto precedente)
#
# 4) Si calcola l'accelerazione dalla forza calcolata in x.m0
#    (si noti come dal punto di vista computazionale sarebbe stato
#     più efficiente calcolare una tantum k/m prima del loop,
#     ma per uno script di questo tipo tale accortezza non è cruciale e
#     rende lo script più intellegibile)
#
# 5-7) Calcoliamo finalmente velocità, posizione e accelerazione 
#      in corrispondenza di t[i]
#