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

save.plots  <- TRUE     # TRUE/FALSE per salvare i plot su file 
metodo      <- 2        # 1: semplice
                        # 2: usa acc di alpha intermedia
dt          <- 0.01     # (s)
tmax        <- 100      # (s)
nmax        <- 100000   # massimo numero di step

# pendolo semplice con massa sospesa di raggio R
g <- 9.8    # (m/s^2)
l <- 1.50   # distanza da punto oscillazione a centro sfera (m)
alpha0   <- 150 # angolo iniziale (gradi)
alpha0 <- alpha0*pi/180
t.range <- c(0,6) # range temporale per plot 
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 * 2/3

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

#----- Functions ---------------------------------------------------
#
#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)*sin(alfa)
   
   return(alfa.pp)
}
# function pausa
pausa <- function() {
  cat ("\n >> Guarda il plot e dai enter per continuare\n")
  scan()
}
# dev.off() o pausa()
dev.off.OR.pause <- function() {
  if (save.plots) {
    dev.off()
  } else {
    pausa()
  }
}
#--------------------------------------------------------------------

n           <- 1
t[n]        <- 0
alpha[n]    <- alpha0
alpha.p[n]  <- alpha.p0
alpha.pp[n] <- acc.ang(alpha[n],alpha.p[n])

# da qui comincia il gioco
while ( t[n] <  tmax && n < nmax) {
  n <- n+1
  t[n] <- t[n-1] + dt
  apm  <- alpha.p[n-1] + acc.ang(alpha[n-1], alpha.p[n-1])*dt/2
  if (metodo == 1) { # metodo rozzo
    alpha.p[n] <- alpha.p[n-1] + acc.ang(alpha[n-1], alpha.p[n-1])*dt
    alpha[n]   <- alpha[n-1] + apm*dt
    
  } else if(metodo == 2) { # metodo leggermente piu' raffinato
    am   <- alpha[n-1] + apm*(dt/2)  # usa accelerazione angolare
                                     # calcolata a meta' intervallino
    alpha.p[n] <- alpha.p[n-1] + acc.ang(am, apm)*dt
    alpha[n]   <- alpha[n-1] + apm*dt    
  }
  alpha.pp[n] <- acc.ang(alpha[n],alpha.p[n])
}

# quindi analizziamo i risultati -------------------------

if (save.plots) {
  nome.base <- sprintf("barra_%.3fm_a%.0fd_M%d_dt%.3f",
                       l, alpha0*180/pi, metodo, dt)
  nome.file <- sprintf("%s__t_alpha.png", nome.base)
  png(filename=nome.file, width=1000, height=600)
}
plot(t,alpha*180/pi, ty='l',   xlab='t (s)', ylab="angolo (gradi)")
dev.off.OR.pause()

if (save.plots) {
  nome.base <- sprintf("barra_%.3fm_a%.0fd_M%d_dt%.3f",
                       l, alpha0*180/pi, metodo, dt)
  nome.file <- sprintf("%s__t_alpha_inizio.png", nome.base)
  png(filename=nome.file, width=1000, height=600)
}
plot(t,alpha*180/pi, ty='l', xlim=t.range, 
     xlab='t (s)', ylab="angolo (gradi)")
dev.off.OR.pause()


if (save.plots) {
  nome.file <- sprintf("%s__t_alpha.p.png", nome.base)
  png(filename=nome.file, width=1000, height=600)
}
plot(t,alpha.p*180/pi, ty='l', xlab='t (s)',
     ylab="velocita' angolare (gradi/s)" )
dev.off.OR.pause()

if (save.plots) {
  nome.file <- sprintf("%s__t_alpha.p_inizio.png", nome.base)
  png(filename=nome.file, width=1000, height=600)
}
plot(t,alpha.p*180/pi, ty='l', xlim=t.range, xlab='t (s)',
     ylab="velocita' angolare (gradi/s)" )
dev.off.OR.pause()


# check alpha.p (scommentare per verificare)
#v.alpha <- rep(0,n-1)
#tv      <- rep(0,n-1)
#for (i in 1:(n-1)) {
#  tv[i]      <- ((i-1) + 1/2)*dt
#  v.alpha[i] <- (alpha[i+1] - alpha[i])/dt
#}
#plot(tv,v.alpha*180/pi, ty='l', xlab='t (s)',
#     ylab="velocita' angolare (gradi/s)")
#pausa()
#----------------------------------------

if (save.plots) {
  nome.file <- sprintf("%s__alpha_alpha.p.png", nome.base)
  png(filename=nome.file, width=1000, height=600)
}
plot(alpha*180/pi,alpha.p*180/pi, ty='l',
     xlab="angolo (gradi)", ylab="velocita' angolare (gradi/s)")
dev.off.OR.pause()

if (save.plots) {
  nome.file <- sprintf("%s__alpha_alpha.p_zoom.png", nome.base)
  png(filename=nome.file, width=1000, height=600)
}
plot(alpha[1:2000]*180/pi,alpha.p[1:2000]*180/pi, ty='l',
     xlab="angolo (gradi)", ylab="velocita' angolare (gradi/s)")
dev.off.OR.pause()