#
#  Quarto di cerchio: valutazione numerica dell'area (e di pi)
#

# funzione del quarto di cerchio
# (generico nome fun per riutilizzarlo con altre funzioni)
fun <- function(x) sqrt(1-x^2)  # att!! senza controlli!!

# per un esempio di altra funzione, ben più complicata, vedi in fondo


N = 30   # nr di intervalli
stima.pi = TRUE

# disegniamo il quarto di cerchio di raggio unitario
x <- seq(0, 1, len=201)         # x usato per il disegno
plot(x, fun(x), ylim=c(0, max(fun(x))), ty='l', asp=1, col='blue')
abline(h=0)
abline(v=0)

x <- seq(0, 1, len=(N+1))       # x usato per discretizzazione
dx = 1 / N            

l = 0
A = 0

for (i in 1:length(x)) {
    xi <- x[i]
    points(rep(xi,2) , c(0, fun(xi) ) , ty='l', lty=2)
    if (xi > 0) {
        xc <- xi - dx/2
        points(xc, fun(xc), pch=19)
        points(c(x[i-1], x[i]), c(fun(x[i-1]), fun(x[i])), ty='l', col='red')

        A = A + dx * fun(xc)
        dy = fun(x[i]) - fun(x[i-1])
        l = l + sqrt(dx^2 + dy^2)
    }    
}

cat(sprintf("N = %d\n", N))

if(stima.pi) {
    cat(sprintf("Area = %f  (esatta: %f)\n", A, pi/4))
    cat(sprintf("Arco = %f  (esatto: %f)\n", l, pi/2))
    cat(sprintf("\nStime di pi greco\n"))
    cat(sprintf("Dall'area: %f\n", 4 * A))    
    cat(sprintf("Dall'arco: %f\n", 2 * l))    
} else {
    cat(sprintf("Area = %f\n", A))    
    cat(sprintf("Arco = %f\n", l))
}


# esempio di funzione complicata
# fun <- function(x)  1 / log( 1 + sin(1+x)^4 ) - 1
# comando da aggiungere dopo plot per avere il titolo
# title(" y = 1 / log( 1 + sin(1+x)^4 ) - 1 ")
# (Att: in questo caso è preferibile togliere 'asp=1')