[Keith Vetter] 2003-02-18 - a Monte Carlo simulation of Buffon's Needle experiment. 

In 1777, Georges Louis Leclerc, Comte de Buffon asked and solved the following problem: what is the probability that a needle of length L when dropped onto a wooden floor with evenly spaced cracks, distance D apart, will land on a crack?

The answer, P = 2L/�D, leads to an amusing way to estimate the value of Pi.

One amusing gotcha in doing the simulation that trips up a lot of people is this: to simulate dropping the needle, one random value you need is the angle the needle is dropped at. If you're not careful, you might code this as '''set angle [[expr 2*$PI*rand()]]'''. You're using Pi to estimate Pi! (If you use degrees you're probably just letting the language library do the conversion to radians for you.) A cute way around this problem is to exploit the fact that Pi is transcendental and just use the consecutive values 1,2,3.... 

[KBK] Hmm, another cute way to draw a random angle is:
   * Choose a random number ''x'', -1 <= ''x'' <= 1.
   * Choose another random number ''y'', -1 <= ''y'' <= 1
   * If ''x''**2 + ''y''**2 > 1, go back to the top.
   * (''x'',''y'') is now a random point inside the unit circle.  Return atan2( ''y'',''x'' )
This method doesn't use Pi to estimate Pi, but does choose a random angle, unlike the
method of simply using integers.

----
 ##+###################################################
 #
 # buffon.tcl -- Simulation of Buffon's Needle Problem
 # by Keith Vetter - Feb 9, 2003
 #
 
 set S(unit) 40                                  ;# Length of needle
 set S(num) 500
 set S(tm) 0
 set S(draw) 0
 set S(total) 0
 set S(cross) 0
 set S(cosn) 1                                   ;# We'll incrementally...
 set S(sinn) 0                                   ;# compute sin and cos
 set S(cos1) [expr {cos(1)}]
 set S(sin1) [expr {sin(1)}]
 
 proc DoDisplay {} {
    wm title . "Buffon's Needle"
    pack [frame .ctrl -relief ridge -bd 2 -padx 5 -pady 5] \
        -side right -fill both -ipady 5
    canvas .c -relief raised -bd 2 -height 500 -width 500 -bg \#ffe4c4
    .c xview moveto 0 ; .c yview moveto 0
    pack .c -side top -fill both -expand 1
 
    DoCtrlFrame
    for {set y $::S(tm)} {$y < 4000} {incr y $::S(unit)} {
        .c create line 0 $y 4000 $y -tag grid
    }
    Clear
    update
 
    trace variable ::S(draw) w Tracer
    set ::S(draw) 0
 }
 proc DoCtrlFrame {} {
    scale .s -from 1 -to 5000 -orient h -variable S(num) -label Needles \
        -relief ridge -highlightthickness 0
    
    button .start  -text "Start" -command Start -bd 4
    button .clear -text "Clear"  -command Clear -bd 4
    button .stop  -text "Stop" -command {set S(draw) 0} -bd 4
    .start configure  -font "[font actual [.start cget -font]] -weight bold"
    .clear configure -font [.start cget -font]
    .stop configure  -font [.start cget -font]
 
    frame .pi -bd 2 -relief raised -pady 5
    label .lbl -text "PI Estimation" -bd 2 -relief groove 
    label .ltot -text "Total"
    label .etot -textvariable S(total) -width 6 -bd 2 -relief sunken -bg white \
        -anchor w
    label .lcross -text "Cross"
    label .ecross -textvariable S(cross) -width 6 -bd 2 -relief sunken -bg white \
        -anchor w
    label .lpi -text "PI"
    label .epi -textvariable S(pi) -width 6 -bd 2 -relief sunken -bg white \
        -anchor w
    button .about -text About -command About
 
    grid .s     - -in .ctrl -sticky ew -row 0
    grid .start - -in .ctrl -sticky ew -row 2
    grid .clear - -in .ctrl -sticky ew
    grid .stop  - -in .ctrl -sticky ew -pady 10
    grid .pi    - -in .ctrl -sticky ew -row 11
    grid .lbl - -in .pi -sticky ew -row 0
    grid .ltot .etot -in .pi -sticky ew -row 2
    grid .lcross .ecross -in .pi -sticky ew
    grid .lpi .epi -in .pi -sticky ew
    
    grid rowconfigure .ctrl 1 -minsize 10
    grid rowconfigure .ctrl 10 -minsize 50
    grid rowconfigure .ctrl 50 -weight 1
    grid columnconfigure .ctrl 1 -weight 1
    grid rowconfigure .pi 1 -minsize 10
 
    grid .about - -in .ctrl -row 100 -sticky ew
 }
 
 proc Clear {} {
    .c delete needle
    set ::S(total) 0
    set ::S(cross) 0
    set ::S(pi) ""
 }
 proc NextSin {} {
    global S
 
    set s [expr {$S(sinn) * $S(cos1) + $S(cosn) * $S(sin1)}]
    set c [expr {$S(cosn) * $S(cos1) - $S(sinn) * $S(sin1)}]
    set S(sinn) $s
    set S(cosn) $c
    return [list [expr {$S(unit) * $s}] [expr {$S(unit) * $c}]]
 }
 proc Start {} {
    set ::S(draw) 1
    DropNeedle $::S(num)
 }
 proc About {} {
    set msg "Buffon's Needle\nby Keith Vetter, February 2003\n\n"
    append msg "A Monte Carlo method to estimate the value of pi invented\n"
    append msg "Comte de Buffon in 1777.\n\n"
    append msg "If you drop a 3 inch needle onto a wooden floor with 3 inch\n"
    append msg "slats, the probability that the needle will cross a crack\n"
    append msg "is 2/pi. Unfortunately, it is very poor way to get an\n"
    append msg "accurate value of pi.\n"
    tk_messageBox -title "About Buffon's Needle" -message $msg
 }
 proc Tracer {var1 var2 op} {
    if {$::S(draw) == 0} {                      ;# Stopping drawing
        .stop config -state disabled
        .start config -state normal
        .clear config -state normal
    } else {                                    ;# Starting to draw
        .stop config -state normal
        .start config -state disabled
        .clear config -state disabled
    }
 }
 proc DropNeedle {{cnt 1}} {
    global S
 
    for {set i 0} {$i < $cnt} {incr i} {
        if {$S(draw) == 0} break
 
        .c itemconfig needle -width 1
        set w [winfo width .c]
        set h [winfo height .c]
        
        # Drop needle between lines visible on the screen
        set x [expr {$w * rand()}]
        set maxy [expr {int(($h - $S(tm))/$S(unit)) * $S(unit)}]
        set y [expr {$S(tm) + $maxy * rand()}]
 
        # Get other end point
        foreach {s c} [NextSin] break
        set x1 [expr {$x + $c}]
        set y1 [expr {$y + $s}]
 
        set yy [expr {($y - $S(tm))}]
        set yy [expr {$yy - $S(unit) * int($yy / $S(unit))}]
        set yy [expr {$yy + $s}]
        
        set color black
        incr S(total)
        if {$yy <= 0 || $yy >= $S(unit)} {      ;# Did it cross
            set color red
            incr S(cross)
        }
        set S(pi) [expr {$S(cross) / .2 / $S(total)}]
        .c create line $x $y $x1 $y1 -fill $color -tag needle -width 5
        if {$i < 50 || ($S(total) % 50) == 0} update
    }
    set S(draw) 0
 }
 DoDisplay

----

[Category Graphics] | [Category Mathematics] | [Category Application] | [Category Whizzlet]