Galton box

A simple implementation of a Galton Box. Balls drop on a lattice of pegs and fall one way or the other; they drop out at the end binomially distributed.

The program takes one argument, the number of balls to drop. The default is 100 balls.

This is the simple form of the implementation, written for clarity.
#
# Draw a galton box, or whatever it's called.
#
# Shamim Mohamed, Sept. 29 1998

link graphics

global pegsize, height, width, pegsize2
   
procedure main(args)
   local n, steps
   steps := 10
   pegsize := 10
   pegsize2 := pegsize * 2
   n := integer(args[1]) | 100
   setup_window(steps)
   every 1 to n do galton(steps)
   WDone()
end

procedure setup_window(n)
   # Draw the n levels of pegs

   # Pegboard size is 2n-1 square
   # Expected max value of histogram is (n, n/2)/2^n 
   # ... approximate with something simpler?
   max := n*n/pegsize
   width := (2*n+1)*pegsize
   height := width + n*n/2*pegsize
   Window("size=" || width || "," || height, "fg=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
      }
   }
   # Set up drawing mode to draw the falling balls
   WAttrib("fg=black")
   WAttrib("drawop=reverse")
end

# Do it!
procedure galton(n)
   local xpos, ypos, i, oldx, oldy
   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
      animate(oldx, oldy, xpos, ypos)
      oldx := xpos
      oldy := ypos
   }
   # Now the ball falls to the floor
   animate(xpos, ypos, xpos, ypos + 40)
   animate(xpos, ypos+40, xpos, ypos + 200)

   # Record this ball
   draw_ball(xpos)
end

procedure animate(xfrom, yfrom, xto, yto)
   animate_actual(xfrom, yfrom, xto, yfrom, 4)
   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)
      WDelay(1)
      FillArc(x, y, pegsize, pegsize)
      x +:= xstep
      y +:= ystep
   }
end

procedure draw_ball(x)
   static ballcounts
   initial ballcounts := table(0)

   ballcounts[x] +:= 1
   FillArc(x, height-ballcounts[x]*pegsize, pegsize, pegsize)
end


Now, a slightly more elaborate one: this one draws a histogram in the background, calculates the size of the window exactly, and allows the animation to be toggled. One problem with the above solution is that with a large number of balls, the animation is fun to watch but takes forever. With this version, hitting the spacebar turns off the animation so you can watch the balls collect more quickly. At any time pressing the spacebar will resume the animation of the balls.

It also accepts an option: galton -pegs 7 50 will simulate 50 balls being dropped through a lattice 7 levels deep.

It's also a little better documented.

 
#
# 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



Shamim Mohamed
Last modified: Mon Oct 5 18:59:47 PDT 1998