#--------------------------------------------------------------
# PED 09/10 (D'Agostini)
# bozza di programma per analizzare i dati del volano
#
# questo non e' un vero script, ma una sequenza di istruzioni
# da eseguire da console in modo interattivo
#--------------------------------------------------------------


save.plots <- TRUE    # sostituire con FALSE se si non vogliono salvare i plots


#--------------------------------------------------------------
# leggiamo i dati dal file NObulNOpala.txt
# note: evitare gli spazi nei nomi dei file  
mis1 <- read.table("NObulNOpala.txt", dec=',', header=TRUE)
# (il parametro dec serve perche' nei dati il decimale
#  era indicata dalla virgola anziche' dal punto)
# per help sul comando:
?read.table


#--------------------------------------------------------------
# info su mis1 (scopriamo che e' una variabile del tipo data.frame):
str(mis1)

# sommario (info statistiche sulle variabili):
summary(mis1)

# primo plot
plot(mis1)

# idem, ma co type line:
plot(mis1, ty='l')

# possiamo visualizzare sia linea che punti:
points(mis1)

# modo alternativo di plottare i dati, dando
# esplicitamente i vettori delle ascisse e delle ordinate
plot(mis1$t, mis1$z)

# ovviamente possiamo anche visualizzare i valori:
mis1$t
mis1$z

# questo ci permette di plottare fino ad un certo tempo,
# ad esempio
plot(mis1$t[mis1$t<10], mis1$z[mis1$t<10])
# stimando cosi' che la discesa dura circa 8.5 secondi

# riproviamo:
plot(mis1$t[mis1$t<8.5], mis1$z[mis1$t<8.5])
# ok per il momento, poi si cerchera' di fare di meglio

# e, ovviamente, si poteva cercare il tempo di inversione stampando
# 'brutalmente' il data frame:
mis1

# per visualizzare i primi 100 valori
mis1[1:100,]

# selezioniamo i dati relativi alla prima discesa e li mettiamo
# nei vettori t e z (diversi da mis1$t !)
t <- mis1$t[mis1$t<8.5]
z <- mis1$z[mis1$t<8.5]

# ( si potevano anche mettere le info selezionate in un
#   altro data frame, mediante
#   sel <- mis1[mis1$t<8,]
# )

plot(t, z)

# plot linearizzato
plot(t^2, z)

# si notano i punti iniziali in cui il sistema e' ancora fermo
# zoomiamo:
plot(t[t<2], z[t<2])

# sembra che il sistema comincia a muoversi con regolarita'
# dal 13.mo punto: selezioniamo i punti da qui in poi:
t <- t[13:length(t)]
z <- z[13:length(z)]
# e riplottiamo
plot(t, z)

# ridefiniamo la scala temporale, prendendo come inizio
# il primo tempo:
t <- t - t[1]
# e riplottiamo
plot(t, z)
# nota: all'inizio l'andamento non e' lineare per effetto
# della velocita' iniziale diversa da zero
# -> come la si puo' valutare?

#------------------------------------------------------------
# Ora passiamo alla velocita'
# prepariamo un vettore vuoto per i valori di velocita'
v <- rep(0, length(t)-1)

# riempiamolo:
for (i in 1:length(v)) {
  v[i] <- (z[i+1] - z[i]) / (t[i+1] - t[i])
}

# plottiamolo (non in funzione del tempo, ma del nr d'ordine):
plot(v)

# ricordiamoci che le velocita' sono calcolate come valori medi
# e quindi si riferiscono ai tempi intermedi:
# calcoliamo i tempi intermedi:
tv <- rep(0, length(v))
for (i in 1:length(tv)) {
  tv[i] <- t[i] + (t[i+1] - t[i])/2
}

# per verificare, visualizziamo t e tv:
tv
t

# ora possiamo plottare la velocita' in funzione del tempo:
plot(tv, v)

# perche' i punti fluttuano cosi' tanto?

# proviamo a far passare qualche retta per i punti
# provando e riprovando (poi cercheremo di fare meglio)
abline(0, -0.35/7)
abline(-0.03, -0.32/7)
abline(-0.04, -0.33/7)
abline(-0.04, -0.32/7)

#--------------------------------------------------------------
#altro modo di calcolare la velocita', usando coppie
#consecutive di punti sperimentali

# nr di coppie (scartando eventuale punto finale spaiato)
nc <- floor(length(z)/2)

vc  <- rep(0, nc)   # nuovi vettori velocita' e tempi
tvc <- rep(0, nc)

for (i in 1:nc) {
  vc[i]  <- (z[2*i] - z[2*i-1]) / (t[2*i] - t[2*i-1])
  tvc[i] <- t[2*i-1] + (t[2*i] - t[2*i-1])/2
}

plot(tvc, vc)

# meglio v o vc?
#--------------------------------------------------------------


#---------------------------------------------------------------
#---------------------------------------------------------------
# ora passiamo all'accelerazione:
# => provare  


#-----------------------------------------------------
# passi successivi:
# 
#  passare da questi comandi 'esplorativi'
#  ad un vero programma di analisi
#  => provare (e riprovare!)
#

#----------------------------------------------------------------
# analisi completa di velocita' e accelerazione
#----------------------------------------------------------------

# velocita' -----------------------------------------------------
nv <- floor(length(mis1[,1])/2)
vc.a  <- rep(0, nv)   # nuovi vettori velocita' e tempi
tvc.a <- rep(0, nv)

for (i in 1:nv) {
  vc.a[i]  <- (mis1$z[2*i] - mis1$z[2*i-1]) / (mis1$t[2*i] - mis1$t[2*i-1])
  tvc.a[i] <- mis1$t[2*i-1] + (mis1$t[2*i] - mis1$t[2*i-1])/2
}

if (save.plots) {
  png(filename="v_tot.png", width=1200, height=800)
}
plot(tvc.a, vc.a, ty='l')
points(tvc.a, vc.a)
if (save.plots) {
  dev.off()
} else {
  cat ("\n >> Guarda il plot e dai enter per continuare\n")
  scan()
}

# accelerazione --------------------------------------------------
na <- floor(length(vc.a)/2)
ac.a  <- rep(0, na)   # nuovi vettori velocita' e tempi
tac.a <- rep(0, na)

for (i in 1:na) {
  ac.a[i]  <- (vc.a[2*i] - vc.a[2*i-1]) / (tvc.a[2*i] - tvc.a[2*i-1])
  tac.a[i] <- tvc.a[2*i-1] + (tvc.a[2*i] - tvc.a[2*i-1])/2
}

if (save.plots) {
  png(filename="a_tot.png", width=1200, height=800)
}
plot(tac.a, ac.a, ty='l')
points(tac.a, ac.a)
grid()
if (save.plots) {
  dev.off()
} else {
  cat ("\n >> Guarda il plot e dai enter per continuare\n")
  scan()
}

#--------- riplottiamo pure le posizioni:
if (save.plots) {
  png(filename="z_tot.png", width=1200, height=800)
}
plot(mis1, ty='l')
points(mis1)
grid()
if (save.plots) {
  dev.off()
} else {
  cat ("\n >> Guarda il plot e dai enter per continuare\n")
  scan()
}

#--------- e le velocita' nella prima discesa:
if (save.plots) {
  png(filename="v_disc1.png", width=1200, height=800)
}
plot(tvc, vc) 
if (save.plots) {
  dev.off()
} else {
  cat ("\n >> Guarda il plot e dai enter per continuare\n")
  scan()
}


# suggerimento fare una function, ad esempio "chiudi.o.pausa()"
# che sostituisce queste 5 righe che si ripetono dopo ogni
# plot