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

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

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

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