# fit lineare fatto a mano
m.molla <- 0.063
m.disco <- 0.0789
n.min   <- 3      # numero minimo di dischi affinche' molla oscilli


# leggiamo tabella
molla <- read.table("molla.dat")
l <- molla$V1
t <- molla$V2

# togliamo gli zeri (richiediamo che entrambi siano > 0:
# non strettamente necessario, ma eserzizio di R) 
ok <- l > 0 & t > 0
l <- l[ok]
t <- t[ok]

np <- length(l)
m  <- n.min:(n.min-1+np)*m.disco + m.molla

#--- funzioni preliminari
my.var <- function (x) mean(x^2) - mean(x)^2
my.sd  <- function(x) sqrt(my.var(x))
my.cov <- function(x,y) mean(x*y) - mean(x)*mean(y)
my.cor <- function(x,y) my.cov(x,y)/(my.sd(x)*my.sd(y)) 
my.skew <- function (x) mean((x-mean(x))^3)/my.sd(x)^3
my.kurt <- function (x) mean((x-mean(x))^4)/my.sd(x)^4
#--- fit lineare
my.lin <- function (x,y) {
  z <- NULL
  z$ok <- FALSE
  if ( length(x) < 2 || length(y) < 2 )       return(z) 
  if ( length(x) != length(y) )               return(z) 
  if ( any(!is.real(x)) || any(!is.real(y)) ) return(z) 

  z$n <- length(x)
  z$Em <- my.cov(x,y)/my.var(x)   # val att slope
  z$Ec <- mean(y) - z$Em*mean(x)  # val att intercetta

  z$r  <- y - z$Ec - z$Em * x     # residui
  if (np > 2) { z$sy <- sqrt(sum(z$r^2)/(z$n-2))  # sigma_y dai residui
              } else { z$sy <- 0 }
  z$Sm <- z$sy / sqrt(z$n) / sqrt(my.var(x)) # sigma(m)
  z$Sc <- z$Sm * sqrt(mean(x^2))       # sigma(c)
  z$rho<- - mean(x) / sqrt(mean(x^2))  # rho(m,c)

  z$ok <  - TRUE
  z
} 

# usiamo my.lin per i due fit
fit.ml <- my.lin(m,l)
fit.mt <- my.lin(sqrt(m),t)

# plot risultati
plot(m,l)
abline(fit.ml$Ec, fit.ml$Em)

plot(sqrt(m),t)
abline(fit.mt$Ec, fit.mt$Em)

plot(m,t)
curve(fit.mt$Ec+fit.mt$Em*sqrt(x),add=TRUE)


# valore atteso di k e g (`stima di...')
Ek <- (2*pi/fit.mt$Em)^2   # k molla
Eg <- fit.ml$Em*Ek         # g gravita'

# incertezze standard
Sk <- Ek * 2 * fit.mt$Sm/fit.mt$Em  
Sg <- Eg * sqrt( 4* (fit.mt$Sm/fit.mt$Em)^2 + (fit.ml$Sm/fit.ml$Em)^2 )

# covarianza e correlazione
cov.kg <- 2^6 * pi^4 * fit.ml$Em / fit.mt$Em^6 * fit.mt$Sm^2
rho.kg <- cov.kg / (Sk*Sg)

# inferenza su k ricondizionata da g=9.8
Ek.g <- Ek + rho.kg * Sk/Sg * (9.8 - Eg)
Sk.g <- Sk * sqrt(1 - rho.kg^2)