#-----------------------------------------------------
# Integrale fatto con MC, con grafica
# (limitatamente al primo lancio)
#
# GdA, aprile 2023
#-----------------------------------------------------

dx    <- c(46, 120, 167)/100       # cm -> m
alpha <- c(11.531, 1.6944, 0.8749) # m^-1

dsdx <- function(x, alpha) sqrt(1 + alpha^2 * x^2)


xi <- seq(0, dx[1], len=101)
plot(xi, dsdx(xi,alpha[1]),  ty='l', col='blue', lwd=2, 
      xaxs = "i", yaxs = "i", ylim=c(0,6), main='Lancio nr. 1', 
      xlab="x  (m)", ylab="ds/dx")

N = 10000
colore <- c('gray', 'red')
cat(sprintf("\nIntegrale fatto con MonteCarlo (%d punti):\n", N))
xr <- runif(N, 0,dx[1])
yr <- runif(N, 0, dsdx(dx[1], alpha[1]))
sotto <- yr <= dsdx(xr, alpha[1])
points(xr,yr, pch=19, cex=0.3, col=colore[sotto+1]) 
integrale <- sum(sotto)/N * dx[1]* dsdx(dx[1], alpha[1]) 
cat(sprintf(" - Lancio %d: %.3f m \n", 1, integrale))