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