log.p = TRUE
curva = TRUE

md = 0.079   # Kg
nd <- 3:10  
x  <- c(3.3, 4.9, 6.8, 8.5, 10.0, 11.8, 13.4, 15.1)
p1 <- c(0.481, 0.539, 0.603, 0.664, 0.692, 0.757, 0.793, 0.846)
p2 <- c(0.449, 0.566, 0.623, 0.693, 0.721, 0.777, 0.806, 0.875)  # att!
p3 <- c(0.491, 0.560, 0.617, 0.672, 0.717, 0.770, 0.810, 0.846)

n <- length(nd)
pm <- (p1+ p3)/3   # senza p2

if(log.p) {
  plot(nd, pm, log='xy', col='blue')
} else {
  plot(nd, pm, col='blue')  
}

# USANDO primo e ultimo punto (NON è la cosa ottimale da fare!)
alpha = ( log(pm[n])-log(pm[1]) ) / (( log(nd[n])-log(nd[1]) ))
print(alpha)
k = pm[1]/nd[1]^alpha
print(k)

# k <- k * 0.995   # aggiustamento a mano (arbitrario!)
# print(k)

if(curva) curve(k*n^alpha, 3, 10, xname='n', col='red', add=TRUE)