 |
#
# Simulate a galton box. The thing that has a network of pegs, and balls
# are dropped in from the top; as it hits each peg it has a 0.5 chance
# of going in either direction. At the bottom, the balls will be
# binomially distributed.
#
# Argument: the number of balls to simulate
# Options: -pegs [# of levels of pegs]
#
# Shamim Mohamed, Sept. 29 1998
link graphics
########################################################################
# globals
global height, width
global pegsize, pegsize2
# pegsize the size of the balls, and also of the pegs
# pegsize2 pre-computed value of pegsize*2 for efficiency
# width, height the dimensions of the window
########################################################################
procedure main(args)
local n, steps
# initialise Icon's random number generator with the time and date
&random := map("sSmMhH", "Hh:Mm:Ss", &clock) +
map("YyXxMmDd", "YyXx/Mm/Dd", &date)
# initial values
steps := 10 # no. of levels in the pegboard
pegsize := 10 # geometric size of the pegs and balls
if *args >= 2 & args[1] == "-pegs" then {
steps := integer(args[2]) | stop("Invalid argument to -pegs")
pop(args); pop(args)
}
n := integer(args[1] | 100) |
stop("Invalid argument for number of balls.")
pegsize2 := pegsize * 2
setup_window(steps, n) # create window; draw the triangle of pegs
every 1 to n do galton(steps)
WDone() # wait for user to type 'q' to exit
end
########################################################################
# Setup the window: draw the n levels of pegs; calculate the expected values
# of the binomial distribution and use that to calculate the window size;
# then draw the binomial distribution as a histogram at the bottom of the
# window.
procedure setup_window(n, nballs)
# Pegboard size is (2n-1)^2
local xpos, ypos, bar_heights, coeff
width := (2*n+1)*pegsize
height := width # placeholder
Window("size=" || width || "," || height, "bg=grayish-white")
WAttrib("fg=dark-grey")
every i := 1 to n do {
ypos := i * pegsize2
xpos := width/2 - (i - 1) * pegsize - pegsize/2
every j := 1 to i do {
FillArc(xpos, ypos, pegsize, pegsize)
xpos +:= pegsize2
}
}
# Calculate the binomial distribution
coeff := 1
h_unit := 1.0 * pegsize * nballs / (2^n)
bar_heights := list(1+(n+1)/2)
every i := 0 to (n+1)/2 do {
# coeff is (n,i)
bar_heights[i+1] := integer(0.5 + coeff * h_unit)
coeff *:= n - i
coeff /:= i + 1
}
# Now we know high big to make the window
height := width + bar_heights[-1] + pegsize2
# this last pegsize*2 being a fudge factor for spacing
WAttrib("size="||width||","||height)
# Draw the binomial bars
xpos := 0
every i := 0 to (n+1)/2 do {
# draw rectangles for (n,i) and (n,n-i)
draw_bar(xpos, bar_heights[i+1])
draw_bar(width - xpos - pegsize, bar_heights[i+1])
xpos +:= pegsize2
}
# Set up drawing mode to draw the falling balls
WAttrib("fg=black")
WAttrib("drawop=reverse")
end
procedure draw_bar(xpos, bar_height)
xpos -:= pegsize/2
DrawLine(xpos, height - bar_height, xpos + pegsize2, height - bar_height)
DrawLine(xpos, height - bar_height, xpos, height)
DrawLine(xpos + pegsize2, height - bar_height, xpos + pegsize2, height)
end
########################################################################
# Do it!
procedure galton(n)
local xpos, ypos, i, oldx, oldy
static quick # animation state; static so it persists
xpos := oldx := width/2 - pegsize/2
ypos := oldy := pegsize
# For every ball...
every 1 to n do {
if ?2 = 1 then
xpos -:= pegsize
else
xpos +:= pegsize
ypos +:= pegsize2
# If spacebar is pressed, toggle animation state
if Event(1) == " " then
if /quick then quick := 1 else quick := &null
if /quick then animate(oldx, oldy, xpos, ypos)
oldx := xpos
oldy := ypos
}
# Now the ball falls to the floor
if /quick then {
animate(xpos, ypos, xpos, ypos + 40)
animate(xpos, ypos+40, xpos, ypos + 200)
}
# Record this ball (it's drawn regardless of animation state)
draw_ball(xpos)
end
procedure animate(xfrom, yfrom, xto, yto)
# We animate in two stages to simulate the ball bouncing off the peg
# It's just a jump to the left
animate_actual(xfrom, yfrom, xto, yfrom, 4)
# And then a step to the-- bottom
animate_actual(xto, yfrom, xto, yto, 10)
end
# Drawing op is already set to "reverse", and fg colour is black.
procedure animate_actual(xfrom, yfrom, xto, yto, steps)
x := xfrom
y := yfrom
xstep := (xto - xfrom)/steps
ystep := (yto - yfrom)/steps
every i := 1 to steps do {
lastx := x
lasty := y
FillArc(x, y, pegsize, pegsize)
delay(1)
FillArc(x, y, pegsize, pegsize)
x +:= xstep
y +:= ystep
}
end
# Fill in one of the buckets at the bottom with the new ball
procedure draw_ball(x)
static ballcounts
initial ballcounts := table(0)
ballcounts[x] +:= 1
FillArc(x, height-ballcounts[x]*pegsize, pegsize, pegsize)
end
|
