#---------------------------------------------------
# Oscillazione della molla: valutazione numerica
#
#  (Per le note a fine linea si veda in fondo)
#
# 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(t, x/max(x), pch=19, col='blue', cex=0.5, xaxs = "i", yaxs = "i",
     xlab='t (s)', ylab='x/xMax;  v/vMax;  a/aMax',
     ylim=c(-1,1.2) )
abline(h=0, col='gray')
T=2*pi/sqrt(k/m)
abline(v=T)
abline(v=2*T)
points(t, v/max(v),  pch=19, col='cyan', cex=0.5)
points(t, a/max(a),  pch=19, col='red', cex=0.5)
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)

#-- 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]
#