# Simuliere inhomogenen Poisson-Prozess radial. 
# Ausdünnung aus Theorem 2.10 

N <- 1000
statistics <- c()

lambda = 2

while(length(statistics) < N){
  cnt <- 0
  R <- 0
  M <- matrix(nrow=2, ncol=2)
  
  while(cnt < 2){
    # Simuliere Punkte des homogenen Poisson-Prozesses
    U <- runif(3)
    R <- R - log(U[1])
    r <- sqrt(R/(pi * lambda))
    
    # Ausdünnung mit Theorem 2.10 
    # (Zur Berechnung der Intensitätsfunktion ist nur der Abstand r des Punktes
    # zum Ursprung relevant)
    
    if(lambda * U[2] <= 1 + exp(-r*r/2)){
      M[1, cnt + 1] <- r * cos(U[3] * 2 * pi)
      M[2, cnt + 1] <- r * sin(U[3] * 2 * pi)
      cnt <- cnt + 1
    }
  }
  
  statistics <- c(statistics, abs(det(M))/2)
  
}

print(mean(statistics))
print(var(statistics))