#
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
# 
#-----

set.seed(0)

# Exercise 32:
#
ellipse_area = function(nSims=10,kappa=1){
  numIn = 0
  x = runif( nSims, min=-1, max=+1 )
  y = runif( nSims, min=-1/kappa, max=+1/kappa )
  ( sum( x^2 + kappa^2 * y^2 < 1 )/nSims ) * ( 4 / kappa ) # return the area 
}

kappas = c(1,0.9,0.75,0.5,0.25,0.1)
ns =  c(1000,5000,10000,50000,1000000)

data = c()
for( k in kappas ){
  for( n in ns ){
    ea = ellipse_area( nSims=n, kappa=k )
    data = c( data, abs(ea - pi/k)/(pi/k) ) # store the error 
  }
}
R = matrix( data, nrow=length(kappas), ncol=length(ns), byrow=TRUE )

rownames(R) = kappas
colnames(R) = ns
print(R)