#-------------------------------------------
# like waiting_2_counts.R, but using R's dgamma()
#  -> note the names of the arguments!   
#
# [Moreover, in this version,
#  - the values have also correct physical dimensions;
#  - the axes start from the origine ]
#
# GdA, Feb 2025
#-------------------------------------------

pausa  <- function() { cat ("\n >> press Enter to continue\n"); scan() }

ft <- function(t,r) r^2*t*exp(-r*t)

r = 1   # 1/s (inverso di secondi!) 
t <- seq(0, 8, len=101)
# curve for k=2, using our function ft() 
plot(t, ft(t,r), ty='l', col='cyan', lwd=1.8, ylim=c(0,1),
     xaxs="i", yaxs="i", xlab= "t (s)", ylab="f(t)")
pausa()
# curve for k=2, using dgamma() 
k = 2
points(t, dgamma(t, shape=k, rate=r),  ty='l', col='blue',  lwd=1.8)
pausa()
# curve for k=1, using dgamma
k = 1
points(t, dgamma(t, shape=k, rate=r),  ty='l', col='red', lty=2, lwd=1.8)

# integrals using our function
cat("\nvalori ottenuti da integrali numerici con la nostra funzione:\n")
E.t <- integrate(function(t) t*ft(t,r), 0, 10)$value
cat(sprintf("  E(t) = %.4f s\n", E.t))
E.t2 <- integrate(function(t) t^2*ft(t,r), 0, 10)$value
cat(sprintf("  sigma(t) = %.4f s\n", sqrt(E.t2-E.t^2)))

# integrals using dgamma
cat("\nvalori ottenuti da integrali numerici con dgamma():\n")
E.t <- integrate(function(t) t*dgamma(t,  shape=2, rate=1), 0, 10)$value
cat(sprintf("  E(t) = %.4f s \n", E.t))
E.t2 <- integrate(function(t) t^2*dgamma(t,  shape=2, rate=1), 0, 10)$value
cat(sprintf("  sigma(t) = %.4f s\n", sqrt(E.t2-E.t^2)))

# Exact values
cat("\nvalori esatti (-> vedi slides):\n")
k=2
E.t <- k/r
cat(sprintf("  E(t) = %.4f s\n", E.t))
Var.t <- k/r^2
cat(sprintf("  sigma(t) = %.4f s\n", sqrt(Var.t)))