pendulum.tcl
1 # pendulum.tcl -- 2 # 3 # This demonstration illustrates how Tcl/Tk can be used to construct 4 # simulations of physical systems. 5 6 if {![info exists widgetDemo]} { 7 error "This script should be run from the \"widget\" demo." 8 } 9 10 package require Tk 11 12 set w .pendulum 13 catch {destroy $w} 14 toplevel $w 15 wm title $w "Pendulum Animation Demonstration" 16 wm iconname $w "pendulum" 17 positionWindow $w 18 19 label $w.msg -font $font -wraplength 4i -justify left -text "This demonstration shows how Tcl/Tk can be used to carry out animations that are linked to simulations of physical systems. In the left canvas is a graphical representation of the physical system itself, a simple pendulum, and in the right canvas is a graph of the phase space of the system, which is a plot of the angle (relative to the vertical) against the angular velocity. The pendulum bob may be repositioned by clicking and dragging anywhere on the left canvas." 20 pack $w.msg 21 22 ## See Code / Dismiss buttons 23 set btns [addSeeDismiss $w.buttons $w] 24 pack $btns -side bottom -fill x 25 26 # Create some structural widgets 27 pack [panedwindow $w.p] -fill both -expand 1 28 $w.p add [labelframe $w.p.l1 -text "Pendulum Simulation"] 29 $w.p add [labelframe $w.p.l2 -text "Phase Space"] 30 31 # Create the canvas containing the graphical representation of the 32 # simulated system. 33 canvas $w.c -width 320 -height 200 -background white -bd 2 -relief sunken 34 $w.c create text 5 5 -anchor nw -text "Click to Adjust Bob Start Position" 35 # Coordinates of these items don't matter; they will be set properly below 36 $w.c create line 0 25 320 25 -tags plate -fill grey50 -width 2 37 $w.c create oval 155 20 165 30 -tags pivot -fill grey50 -outline {} 38 $w.c create line 1 1 1 1 -tags rod -fill black -width 3 39 $w.c create oval 1 1 2 2 -tags bob -fill yellow -outline black 40 pack $w.c -in $w.p.l1 -fill both -expand true 41 42 # Create the canvas containing the phase space graph; this consists of 43 # a line that gets gradually paler as it ages, which is an extremely 44 # effective visual trick. 45 canvas $w.k -width 320 -height 200 -background white -bd 2 -relief sunken 46 $w.k create line 160 200 160 0 -fill grey75 -arrow last -tags y_axis 47 $w.k create line 0 100 320 100 -fill grey75 -arrow last -tags x_axis 48 for {set i 90} {$i>=0} {incr i -10} { 49 # Coordinates of these items don't matter; they will be set properly below 50 $w.k create line 0 0 1 1 -smooth true -tags graph$i -fill grey$i 51 } 52 53 $w.k create text 0 0 -anchor ne -text "\u03b8" -tags label_theta 54 $w.k create text 0 0 -anchor ne -text "\u03b4\u03b8" -tags label_dtheta 55 pack $w.k -in $w.p.l2 -fill both -expand true 56 57 # Initialize some variables 58 set points {} 59 set Theta 45.0 60 set dTheta 0.0 61 set pi 3.1415926535897933 62 set length 150 63 set home 160 64 65 # This procedure makes the pendulum appear at the correct place on the 66 # canvas. If the additional arguments "at $x $y" are passed (the 'at' 67 # is really just syntactic sugar) instead of computing the position of 68 # the pendulum from the length of the pendulum rod and its angle, the 69 # length and angle are computed in reverse from the given location 70 # (which is taken to be the centre of the pendulum bob.) 71 proc showPendulum {canvas {at {}} {x {}} {y {}}} { 72 global Theta dTheta pi length home 73 if {$at eq "at" && ($x!=$home || $y!=25)} { 74 set dTheta 0.0 75 set x2 [expr {$x - $home}] 76 set y2 [expr {$y - 25}] 77 set length [expr {hypot($x2, $y2)}] 78 set Theta [expr {atan2($x2, $y2) * 180/$pi}] 79 } else { 80 set angle [expr {$Theta * $pi/180}] 81 set x [expr {$home + $length*sin($angle)}] 82 set y [expr {25 + $length*cos($angle)}] 83 } 84 $canvas coords rod $home 25 $x $y 85 $canvas coords bob \ 86 [expr {$x-15}] [expr {$y-15}] [expr {$x+15}] [expr {$y+15}] 87 } 88 showPendulum $w.c 89 90 # Update the phase-space graph according to the current angle and the 91 # rate at which the angle is changing (the first derivative with 92 # respect to time.) 93 proc showPhase {canvas} { 94 global Theta dTheta points psw psh 95 lappend points [expr {$Theta+$psw}] [expr {-20*$dTheta+$psh}] 96 if {[llength $points] > 100} { 97 set points [lrange $points end-99 end] 98 } 99 for {set i 0} {$i<100} {incr i 10} { 100 set list [lrange $points end-[expr {$i-1}] end-[expr {$i-12}]] 101 if {[llength $list] >= 4} { 102 $canvas coords graph$i $list 103 } 104 } 105 } 106 107 # Set up some bindings on the canvases. Note that when the user 108 # clicks we stop the animation until they release the mouse 109 # button. Also note that both canvases are sensitive to <Configure> 110 # events, which allows them to find out when they have been resized by 111 # the user. 112 bind $w.c <Destroy> { 113 after cancel $animationCallbacks(pendulum) 114 unset animationCallbacks(pendulum) 115 } 116 bind $w.c <Button-1> { 117 after cancel $animationCallbacks(pendulum) 118 showPendulum %W at %x %y 119 } 120 bind $w.c <B1-Motion> { 121 showPendulum %W at %x %y 122 } 123 bind $w.c <ButtonRelease-1> { 124 showPendulum %W at %x %y 125 set animationCallbacks(pendulum) [after 15 repeat [winfo toplevel %W]] 126 } 127 bind $w.c <Configure> { 128 %W coords plate 0 25 %w 25 129 set home [expr {%w/2}] 130 %W coords pivot [expr {$home-5}] 20 [expr {$home+5}] 30 131 } 132 bind $w.k <Configure> { 133 set psh [expr {%h/2}] 134 set psw [expr {%w/2}] 135 %W coords x_axis 2 $psh [expr {%w-2}] $psh 136 %W coords y_axis $psw [expr {%h-2}] $psw 2 137 %W coords label_dtheta [expr {$psw-4}] 6 138 %W coords label_theta [expr {%w-6}] [expr {$psh+4}] 139 } 140 141 # This procedure is the "business" part of the simulation that does 142 # simple numerical integration of the formula for a simple rotational 143 # pendulum. 144 proc recomputeAngle {} { 145 global Theta dTheta pi length 146 set scaling [expr {3000.0/$length/$length}] 147 148 # To estimate the integration accurately, we really need to 149 # compute the end-point of our time-step. But to do *that*, we 150 # need to estimate the integration accurately! So we try this 151 # technique, which is inaccurate, but better than doing it in a 152 # single step. What we really want is bound up in the 153 # differential equation: 154 # .. - sin theta 155 # theta + theta = ----------- 156 # length 157 # But my math skills are not good enough to solve this! 158 159 # first estimate 160 set firstDDTheta [expr {-sin($Theta * $pi/180)*$scaling}] 161 set midDTheta [expr {$dTheta + $firstDDTheta}] 162 set midTheta [expr {$Theta + ($dTheta + $midDTheta)/2}] 163 # second estimate 164 set midDDTheta [expr {-sin($midTheta * $pi/180)*$scaling}] 165 set midDTheta [expr {$dTheta + ($firstDDTheta + $midDDTheta)/2}] 166 set midTheta [expr {$Theta + ($dTheta + $midDTheta)/2}] 167 # Now we do a double-estimate approach for getting the final value 168 # first estimate 169 set midDDTheta [expr {-sin($midTheta * $pi/180)*$scaling}] 170 set lastDTheta [expr {$midDTheta + $midDDTheta}] 171 set lastTheta [expr {$midTheta + ($midDTheta + $lastDTheta)/2}] 172 # second estimate 173 set lastDDTheta [expr {-sin($lastTheta * $pi/180)*$scaling}] 174 set lastDTheta [expr {$midDTheta + ($midDDTheta + $lastDDTheta)/2}] 175 set lastTheta [expr {$midTheta + ($midDTheta + $lastDTheta)/2}] 176 # Now put the values back in our globals 177 set dTheta $lastDTheta 178 set Theta $lastTheta 179 } 180 181 # This method ties together the simulation engine and the graphical 182 # display code that visualizes it. 183 proc repeat w { 184 global animationCallbacks 185 186 # Simulate 187 recomputeAngle 188 189 # Update the display 190 showPendulum $w.c 191 showPhase $w.k 192 193 # Reschedule ourselves 194 set animationCallbacks(pendulum) [after 15 [list repeat $w]] 195 } 196 # Start the simulation after a short pause 197 set animationCallbacks(pendulum) [after 500 [list repeat $w]]