# 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) | 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 == "-pegs" then { steps := integer(args) | stop("Invalid argument to -pegs") pop(args); pop(args) } n := integer(args | 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