# R Program to simulate Buffon's Needle

#   by Jeffrey Rosenthal, www.probability.ca

# Source this file in R to set up and drop the first 5 needles.
# Then, type buffon(n) to add n more needles (default=5).
# Or, type buffon(n,add=FALSE) to start over with n fresh needles.


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

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

# Draw the background.
if ( (!add) | is.null(x.global) ) {
  x.start = y.start = angle = NULL
  plot(0, 0, xlim=c(0,size), ylim=c(0,size), type='n', xlab="", ylab="")
  segments(rep(-1,size+1), 0:size, rep(size+1,size+1), 0:size, col="blue")
} else {
  x.start = x.global
  y.start = y.global
  angle = an.global
}

# Generate the new needle min-points and angles.
firstnew = length(x.start) + 1
x.start = c(x.start, runif(n, -0.5, size))
y.start = c(y.start, runif(n, 0, size))
angle = c(angle, runif(n, 0, 3.1415926))
n = length(x.start)
x.global <<- x.start
y.global <<- y.start
an.global <<- angle

# Compute the needle endpoints.
x.end = x.start + cos(angle)
y.end = y.start + sin(angle)

# Draw the needles.
segments(x.start[firstnew:n], y.start[firstnew:n],
	x.end[firstnew:n], y.end[firstnew:n], col="red")

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

# Calculate and output 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")

}

# Run it once, to get started.
buffon()