# bozza script per la soluzione numerica
# del moto del pendolo
#
# G. D'Agostini, ottobre 09 (PED 2009/2010)
#
#---------------------------------------------

# pendolo semplice con massa sospesa di raggio R
g <- 9.8    # (m/s^2)
l <- 1.00   # distanza da punto oscillazione a centro sfera (m)
R <- 0.02   # raggio sfera (m)
m <- 0.30   # massa sfera (kg)
beta <- 0.  # coefficiente attrito viscosita' (N/(m/s))
            # per ora lo ignoriamo
alpha0   <- 5 # angolo iniziale (gradi)
alpha0 <- alpha0*pi/180
alpha.p0 <- 0 # velocita' angolare iniziale (gradi/s)
              # [potremmo pensare che la sferetta viene inizialmente
              #  spinta invece che semplicemente lasciata andare]
alpha.p0 <- alpha.p0*pi/180

# lunghezza equivalente per pendolo di massa puntiforme
l.eq <- l * (1 + 2/5 * (R/l)^2)

# arrays in cui conserviamo grandezze di interesse
# ('p' e 'pp' stanno per 'punto' e '2 punti'):
t        <- numeric()   # tempo
alpha    <- numeric()   # angolo (radianti!)
alpha.p  <- numeric()   # vel. angolare (rad/s)
alpha.pp <- numeric()   # acc. angolare (rad/s^2)

#funzione che ci da' l'accelerazione angolare:
acc.ang <- function(alfa, alfa.p) {
   # alfa e alfa.p (in rad e rad/s) sono le variabili locali
   # usate all'interno della function

   # nota: le costanti g, m, R e beta non vengono
   # passate in quanto le function di 'R' hanno accesso
   # alle variabili globali

   # caso elementare, con approssimazione di piccoli angoli
   alfa.pp <- - (g/l.eq)*alfa
   
   return(alfa.pp)
}

dt          <- 0.01   # (s)
nmax        <- 100000 # massimo numero di step
n           <- 1
alpha[n]    <- alpha0
alpha.p[n]  <- alpha.p0
alpha.pp[n] <- acc.ang(alpha[n],alpha.p[n])

# da qui comincia il gioco: completare
while ( .... ) {
  n <- n+1
  ........
}

# quindi analizzare i risultati
plot(...)