#-----------------------------------------------------------
# Soluzione numerica del pendolo usando l'energia potenziale
# (limitatamente a un quarto di periodo)
#
# GdA, aprile 2024
#------------------------------------------------------------

pausa <- function() { cat ("\n >> press Enter to continue\n"); scan() }

# parametri del pendolo  ---------------------------------------------------
l     <- 1                    # m   
g     <- 9.8                  # m/s^2
th0.g <- 5                    # gradi
g2r   <- pi/180               # conversione gradi -> radianti
th0   <- th0.g * g2r          # -> radianti

# pulsazione, periodo e quarto di periodo per 'piccole oscillazioni'
omega0 <- sqrt(g/l)
T0     <- 2*pi/omega0
T0.q   <- T0/4 

# metodo numerico, a step ---------------------------------------------------
n    <- 100                       # numero di intervallini 
Dth  <- th0/n                     # intervallini in theta  
th   <- rev(seq(0, th0, by=Dth))  # valori di theta, da th0 a 0
Ds   <- l*Dth                     # intervallini in s, coordinata curvilinea
h    <- l*(1-cos(th))             # altezze ai vari theta
v    <- sqrt( 2*g*(h[1]-h) )      # v ad ogni altezza, usando conservazione energia
t    <- rep(0, n+1)               # mettiamo tutti zeri in 't'
                                  # [se gli intervallini sono 'n' i valori sono 'n+1']
# loop su tutti i valori 
for (i in 2:(n+1)) {              # Att: in R l'indice dei vettori comincia da 1          
   vm   <- (v[i-1]+v[i])/2        # velocità media nell'intervallino 
   Dt   <- Ds / vm                # tempo per percorrere l'intervallino
   t[i] <- t[i-1] + Dt            # tempo trascorso fino al 'i'
}
T.q  <- t[n+1]                    # quarto di periodo ottenuto numericamente

# risultati -----------------------------------------------------------------
cat(sprintf(" angolo massimo (gradi):                     %.1f\n", th0.g))
cat(sprintf(" numero di intervallini:                     %d\n", n))
cat(sprintf(" T/4 (numerico):                             %.6f\n", T.q))
cat(sprintf(" T/4 (approssimazione per 'piccoli angoli'): %.6f\n", T0.q))
cat(sprintf(" rapporto (numerico su approssimato):        %.6f (=1/%.6f)\n",
            T.q/T0.q, T0.q/T.q))

# plot di theta vs t  (theta in gradi)
plot(t, th/g2r, xlab='t (s)', ylab='theta (gradi)', col='blue')
grid()
# confronto con l'approssimazione per 'piccoli angoli'
points(t, th0.g*cos(omega0*t), ty='l', lty=2, lwd=1.8, col='red')
abline(v=T0.q, col='red', lty=3)
pausa()

# plot di v vs t
plot(t,v,  xlab='t (s)',  ylab='v (m/s)', col='blue',
     ylim=c(0, max(v)*1.1))
grid()
# confronto con l'approssimazione per 'piccoli angoli'
points(t, omega0*th0*sin(omega0*t), ty='l', lty=2, lwd=1.8, col='red')
abline(v=T0.q, col='red', lty=3)
pausa()

# plot di theta vs v   (theta in gradi)
plot(th/g2r, v,  xlab='theta (gradi)', ylab='v (m/s)', col='blue',
     ylim=c(0, max(v)*1.1))
grid()
# confronto con l'approssimazione per 'piccoli angoli'
points(th0.g*cos(omega0*t), omega0*th0*sin(omega0*t), ty='l', lwd=1.8, lty=2, col='red')