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