# buffon() -- a simple R program to simulate the Buffon's Needle experiment

# Inspired by the version at:
# https://www.reddit.com/r/statistics/comments/2abaoe/i_wrote_a_simulation_with_graphics_for_buffons/

# Usage: buffon(), or buffon(100,add=FALSE) etc.

x.gl = y.gl = an.gl = NULL
size = 10

buffon = function(n=10, add=TRUE) {

# Draw the background.
if ( (!add) | (!exists('x.gl')) ) {
  x.mid = y.mid = angle = NULL
  plot(0, 0, xlim=c(0,size), ylim=c(0,size), type='n')
  for (i in 0:size)
    lines (y=c(i,i), x=c(-1,size+1), col="blue")
} else {
  x.mid = x.gl
  y.mid = y.gl
  angle = an.gl
}

# Generate the new needle mid-points and angles.
prevn = length(x.mid)
x.mid = c(x.mid, runif(n,0,size))
y.mid = c(y.mid, runif(n,0,size))
angle = c(angle, runif(n,0,pi/2))
n = length(x.mid)
x.gl <<- x.mid
y.gl <<- y.mid
an.gl <<- angle

# cat(n,x.mid,angle,"\n")

# Compute the needle endpoints.
x.max = x.mid + 0.5* cos(angle)
x.min = x.mid - 0.5* cos(angle)
y.max = y.mid + 0.5* sin(angle)
y.min = y.mid - 0.5* sin(angle)

# draw the needles
for (i in prevn:n)
  lines(x=c(x.min[i],x.max[i]), y=c(y.min[i],y.max[i]), col="red")

# Count the hits.
hits = 0
for (i in 0:size)
  hits = hits + sum( (y.min<i) & (y.max>i) )

# calculate and display the pi estimate
pi.estimate = 2 * n / hits
cat(n, "needles;", hits, "hits;",
	"estimate = 2 x", n, "/", hits, "=", pi.estimate,
	"; error =", pi.estimate-pi, "\n")

}  # End of buffon() function.


# Examples:
buffon(5, add=FALSE)
buffon(50, add=TRUE)
buffon(100)