#-------------------------------------------------
# Propagazione di incertezze mediante Monte Carlo 
# (dati in linearizzazione_propagazioni.pdf)

# GdA maggio 2025
#-------------------------------------------------

# function pausa
pausa <- function() {
  cat ("\n >> Guarda il plot e dai enter per continuare\n")
  scan()
}

# grandezza di 'input'
E.a <- 29.73  # cm    (valore atteso)
s.a <-  0.03  # cm    (incertezza standard)
E.b <- 21.45  # cm
s.b <-  0.04  # cm 

n = 100000
# valori (pseudo-)casuali di a e b, assumendo distr. 'normale'
a <- rnorm(n, E.a, s.a)
b <- rnorm(n, E.b, s.b)

# grandezze derivate 
p <- 2*(a+b)
A <- a*b
d <- sqrt(a^2+b^2)

# plot (istogramma) e summaries 
E.p <- mean(p)
s.p <- sd(p)
ris.p <- sprintf("\n p = %.2f +- %.2f cm\n", E.p, s.p)
cat(ris.p)
h.p <- hist(p, nc=100, col='cyan', xlab='p (cm)')
text(E.p+2.5*s.p, 0.9*max(h.p$counts), cex=2, ris.p)

pausa()
E.A <- mean(A)
s.A <- sd(A)
ris.A <- sprintf("\n A = %.1f +- %.1f cm^2\n", E.A, s.A)
cat(ris.A)
h.A <- hist(A, nc=100, col='cyan', xlab='A (cm^2)')
text(E.A+2.5*s.A, 0.9*max(h.A$counts), cex=2, ris.A)

pausa()

plot(p, A, xlab='p (cm)', ylab='A (cm^2)', col='cyan')
text(E.p - 2.5*s.p, E.A, cex=2, ris.A)
text(E.p, E.A - 2.5*s.A, cex=2, ris.p)
points(E.p, E.A, cex=3)
ris.rho <- sprintf("\n rho(p, A) = %.4f\n", cor(p, A))
cat(ris.rho)
text(E.p, E.A + 3*s.A, cex=2, ris.rho)