# script per lo studio di fattibilita' dell'esperimento
# del piano inclinato per misura di g e di mu_D
# mediante misure di tempi di discesa in funzione dell'angolo
# di inclinazione
#
#  PED 2009/2010  
#                     G. D'Agostini, ottobre 09
#---------------------------------------------------------------

l = 2       # m
g = 9.8     #
mu.d = 0.1  # coeff. attrito dinamico

alpha <- c(10,20,30,40,50,60)*pi/180

#---- function tempo di arrivo  
t.disc <- function(alpha) sqrt (2*l / (g*( sin(alpha) -mu.d*cos(alpha))  ) )

td <- t.disc(alpha)

plot(tan(alpha), l/td^2/cos(alpha) )

# facciamo un fit dell'andamento linearizzato
fit.t <- lm(  l/td^2/cos(alpha) ~  tan(alpha))
fit.t
abline(fit.t, col='red')

# per avere maggiori dettagli sull'oggetto fit.t
str(fit.t)

#----- velocita' finale da tempo di discesa (assumendo
#      moto uniformemente accelerato):
vf <- (l/td)/2

# analisi usando v in funzione di alpha:
plot(tan(alpha), vf^2/cos(alpha) )
fit.v <- lm (vf^2/cos(alpha)~ tan(alpha))
fit.v
abline(fit.v)

#--------------------------------------------------
# torniamo alle misure 
# misure di tempo, aggiungiamo errori di misura:
sigma.t <- 0.05    # sigma piu' o meno 'inventata'
td.e <- td + rnorm(length(td), 0, sigma.t)  # errore gaussiano
plot(tan(alpha), l/td.e^2/cos(alpha) )
fit.t.e <- lm(  l/td.e^2/cos(alpha) ~  tan(alpha))
fit.t.e
abline(fit.t.e, col='red')

# simuliamo qualcosa vicino a quanto fatto in lab., con
# tre misure per ogni angolo
td.e3   <- rep(0, length(alpha)*3)
alpha3  <- rep(0, length(alpha)*3)
for (i in 1:length(alpha) ){
  for (j in 1:3) {
    alpha3[(i-1)*3 + j] <- alpha[i]
    td.e3[(i-1)*3 + j]  <- td[i] + rnorm(1, 0, sigma.t)
  }
}

plot(tan(alpha3), l/td.e3^2/cos(alpha3) )
fit.t.e3 <- lm(  l/td.e3^2/cos(alpha3) ~  tan(alpha3))
fit.t.e3
abline(fit.t.e3, col='red')

# ricaviamoci g e mu.d dai 'dati sperimentali'
g.mis    <- fit.t.e3$coefficients[2]*2
mu.d.mis <- - fit.t.e3$coefficients[1]/fit.t.e3$coefficients[2]

# per l'analisi dei dati sperimentali e' sufficiente
# rimpiazzare alpha3 e td.e3 con i valori sperimentali.

# nota: quando si hanno dubbi su cosa faccia una funzione, si usi l'help,
# ad es.
?lm


#################################################################
#
# facciamo una stima piu' "seria" dei parametri

#---- function tempo di arrivo: nuova funzione nella quale
#---- anche g e mu.d sono dati come parametri 
t.disc1 <- function(alpha, g, mu.d) sqrt (2*l / (g*( sin(alpha) -mu.d*cos(alpha))  ) )

# somma dei quadrati degli scarti fra
# dati sperimentali e andamento teorico
sq <- function(x) {
   return(sum(
              (td.e3 - t.disc1(alpha3, x[1], x[2]))^2
             )
         )
}


# 'non linear minimization' (vedi help)
# usando come starting point i valori ottenuti rozzamente
# dal 'linear model', previa linearizzazione:
fit.mq <-  nlm(sq, c(g.mis, mu.d.mis))

# parametri che minimizzano sq
g.mq    <- fit.mq$estimate[1]
mu.d.mq <- fit.mq$estimate[2]

# guardiamoli
g.mq
mu.d.mq

# plot di punti sperimentali Versus curva teorica
# con i parametri della minimizzazione:
plot(alpha3, td.e3, main='Simulazione esperimento',xlab='alpha', ylab='t')

alpha.th <- seq(from=alpha3[1], to=alpha3[length(alpha3)], len=100)
points(alpha.th, t.disc1(alpha.th, g.mq, mu.d.mq),ty='l', col='blue')
legend("topright", c("punti sperimentali","fit"),
           lty=c(0,1), pch=c(1,-1), col=c("black","blue"))


# PENSARCI SU ....
#
# Inoltre: questo script rappresenta, come al solito,
# una lista di comandi nati da una sessione interattiva
# -> andrebbe ripulito un po'