Version 21 of TclSpringies : A simple mass and spring simulator

Updated 2004-12-08 11:49:01

if 0 {

MBS After visiting the Runge-Kutta page recently, I decided to re-visit something that I was playing with about a year ago...

I got the idea of trying to create a simple mass and spring simulation based on Xspringies [L1 ] after finding a Java based version [L2 ].

I wasn't interested in all the 'bells and whistles' that were included in the Xspringies. I wanted to keep it 'simple'.

The following is what I ended up with. Comments and code fixes are very much welcome.


escargo 6 Dec 2004 -- On Windows XP Pro, clicking the X in the title bar to close the application fails because "file_quit" is not defined (in proc init_display).

MG 7 Dec 2004 - Added a call to "catch" into the wm protocol to "fix" that. I didn't want to remove the wm protocol completely in case MBS was planning to add it in later/had left it out by mistake.

MBS 7 Dec 2004 - I cut-n-pasted that from another app and never noticed that I was still using "file_quit".

AET 7dec04 - What a great little app to demo Tcl/Tk! Well done. The 700x700 canvas is a little large for the 800x600 monitor on the Win2000 machine I tried it on, so all the control buttons disappeared off the bottom of the screen. Alt-click doesn't work on Windows, so I couldn't move the window enough to see them. I changed the canvas to 700x500 when I pasted it into tclkit, which was fine. You may consider putting the controls on the top if you want to keep the canvas size.

MBS 7 Dec 2004 - I have gotten so used to running at a higher resolution, that I never gave it a thought that 700x700 might be a problem. Let me think about it a little to see if I can come up with a "simple / general" type of fix. In the mean time, if anyone has a suggestion or fix, let me know.

escargo 7 Dec 2004 -- On Windows XP Pro, making the top-level window smaller hid the controls Since these controls included the Quit button, and the X on the title bar doesn't work, there was no way to kill the application. You might want to change the packing order.

SS 8 Dec 2004 -- Comments? Very very cool! That should be added to the Tk demo distribution to show what some line of Tcl can do (in the hands of the smart programmer), and how powerful the canvas widget is. And also that Tcl is *not* too slow even for realtime simulation sometimes.


 }

 if 0 {

Initialize global values...

I don't remember why I called it 'State', but State is a global array that holds almost everything.

 }

 global State

 set State(t)    0.0           ;# current time
 set State(dt)   0.1           ;# time step
 set State(MaxTime) 1000.      ;# Max time value

 # After each time step, the maximum distance that any node travels
 # is calculated.  If this distance, "dmax", is less that "tol", then
 # stop the simulation.

 set State(tol)  0.01

 set State(gravity) -9.8
 set State(curNode) none


 if 0 {

Initialize the geometry to be displayed.

node_def is a list of nodes containing the following for each node:

  • node name
  • X coordinate
  • Y coordinate
  • node mass
  • free/fixed (whether that node is free to move or fixed)
  • spring list (the list of springs attached to this node)

spring_def is the list of springs containing the following info for each spring:

  • spring name
  • spring constant : k
  • spring damping constant : kd
  • rest length for this spring
  • node list (the list of nodes attached to this spring)

init_geom has multiple definitions for node_def and spring_def. Each one is just a different test. Select the one you wish to try. The current setting for node_def and spring_def is sort of a simple suspension bridge type of configuration.

I probably should have added some code so that the user could load node_def and spring_def from a file, but I never got around to doing it...

Oh well.

 }

 proc init_geom { } {
     global State


 ########################################################################

     set node_def { \
                        { n1 250.  350.  50. fixed { s1    } } \
                        { n2 300.  350.  50. free  { s1 s2 } } \
                        { n3 250.  300.  50. free  { s2    } } \
                    }

     set spring_def { \
                        { s1 10.0 50.0  50. {n1 n2} } \
                        { s2 10.0 50.0  50. {n2 n3} } \
                      }

 ########################################################################

     set node_def { \
                        { n1  100 350  50. fixed { s1  s5      } } \
                        { n2  150 350  50. free  { s1  s2  s6  } } \
                        { n3  200 350  50. free  { s2  s3  s7  } } \
                        { n4  250 350  50. free  { s3  s4  s8  } } \
                        { n5  300 350  50. fixed { s4  s9      } } \
                        { n6  100 300  50. fixed { s5  s10     } } \
                        { n7  150 300  50. free  { s6  s10 s11 } } \
                        { n8  200 300  50. free  { s7  s11 s12 } } \
                        { n9  250 300  50. free  { s8  s12 s13 } } \
                        { n10 300 300  50. fixed { s9  s13     } } \
                    }

     set spring_def { \
                          { s1  50.0  45.0  50. {n1   n2 } } \
                          { s2  50.0  45.0  50. {n2   n3 } } \
                          { s3  50.0  45.0  50. {n3   n4 } } \
                          { s4  50.0  45.0  50. {n4   n5 } } \
                          { s5  50.0  45.0  50. {n1   n6 } } \
                          { s6  50.0  45.0  50. {n2   n7 } } \
                          { s7  50.0  45.0  50. {n3   n8 } } \
                          { s8  50.0  45.0  50. {n4   n9 } } \
                          { s9  50.0  45.0  50. {n5   n10} } \
                          { s10 50.0  45.0  50. {n6   n7 } } \
                          { s11 50.0  45.0  50. {n7   n8 } } \
                          { s12 50.0  45.0  50. {n8   n9 } } \
                          { s13 50.0  45.0  50. {n9   n10} } \
                      }

 ########################################################################

     set node_def { \
                        { n1  50 350 100. fixed {      s1 } } \
                        { n2 100 350  50. free  { s1   s2 } } \
                        { n3 150 350 100. free  { s2   s3 } } \
                        { n4 200 350  50. free  { s3   s4 } } \
                        { n5 250 350  60. free  { s4   s5 } } \
                        { n6 350 350  70. free  { s5   s6 } } \
                        { n7 450 350 100. fixed { s6      } } \
                    }

     set spring_def { \
                          { s1  30.0  25.0  25. {n1   n2} } \
                          { s2  40.0  35.0  25. {n2   n3} } \
                          { s3  50.0  45.0  25. {n3   n4} } \
                          { s4  50.0  45.0  25. {n4   n5} } \
                          { s5  40.0  35.0  25. {n5   n6} } \
                          { s6  30.0  25.0  25. {n6   n7} } \
                      }

 ########################################################################

     set node_def { \
                        { n1  100 350  50. fixed { s1  s5  s14         } } \
                        { n2  150 350  50. free  { s1  s2  s6  s15 s16 } } \
                        { n3  200 350  50. free  { s2  s3  s7  s17 s18 } } \
                        { n4  250 350  50. free  { s3  s4  s8  s19 s20 } } \
                        { n5  300 350  50. fixed { s4  s9  s21         } } \
                        { n6  100 300  50. fixed { s5  s10 s15         } } \
                        { n7  150 300  50. free  { s6  s10 s11 s14 s17 } } \
                        { n8  200 300  50. free  { s7  s11 s12 s16 s19 } } \
                        { n9  250 300  50. free  { s8  s12 s13 s18 s21 } } \
                        { n10 300 300  50. fixed { s9  s13 s20         } } \
                    }

     set spring_def { \
                          { s1  90.0  45.0  50. {n1   n2 } } \
                          { s2  90.0  45.0  50. {n2   n3 } } \
                          { s3  90.0  45.0  50. {n3   n4 } } \
                          { s4  90.0  45.0  50. {n4   n5 } } \
                          { s5  90.0  45.0  50. {n1   n6 } } \
                          { s6  90.0  45.0  50. {n2   n7 } } \
                          { s7  90.0  45.0  50. {n3   n8 } } \
                          { s8  90.0  45.0  50. {n4   n9 } } \
                          { s9  90.0  45.0  50. {n5   n10} } \
                          { s10 90.0  45.0  50. {n6   n7 } } \
                          { s11 90.0  45.0  50. {n7   n8 } } \
                          { s12 90.0  45.0  50. {n8   n9 } } \
                          { s13 90.0  45.0  50. {n9   n10} } \
                          { s14 90.0  45.0  50. {n1   n7 } } \
                          { s15 90.0  45.0  50. {n2   n6 } } \
                          { s16 90.0  45.0  50. {n2   n8 } } \
                          { s17 90.0  45.0  50. {n3   n7 } } \
                          { s18 90.0  45.0  50. {n3   n9 } } \
                          { s19 90.0  45.0  50. {n4   n8 } } \
                          { s20 90.0  45.0  50. {n4   n10} } \
                          { s21 90.0  45.0  50. {n5   n9 } } \
                      }

 ########################  bridge  ######################################

     set node_def { \
                        { n1   25 300  20. fixed {     s1  s49 } } \
                        { n2   50 300  20. free  { s1  s2  s50 } } \
                        { n3   75 300  20. free  { s2  s3  s51 } } \
                        { n4  100 300  20. free  { s3  s4  s52 } } \
                        { n5  125 300  20. free  { s4  s5  s53 } } \
                        { n6  150 325  20. free  { s5  s6  s54 } } \
                        { n7  175 350  20. free  { s6  s7  s55 } } \
                        { n8  200 375  20. free  { s7  s8  s56 } } \
                        { n9  225 400  20. fixed { s8  s9  s57 } } \
                        { n10 250 380  20. free  { s9  s10 s58 } } \
                        { n11 275 370  20. free  { s10 s11 s59 } } \
                        { n12 300 360  20. free  { s11 s12 s60 } } \
                        { n13 325 350  20. free  { s12 s13 s61 } } \
                        { n14 350 360  20. free  { s13 s14 s62 } } \
                        { n15 375 370  20. free  { s14 s15 s63 } } \
                        { n16 400 380  20. free  { s15 s16 s64 } } \
                        { n17 425 400  20. fixed { s16 s17 s65 } } \
                        { n18 450 375  20. free  { s17 s18 s66 } } \
                        { n19 475 350  20. free  { s18 s19 s67 } } \
                        { n20 500 325  20. free  { s19 s20 s68 } } \
                        { n21 525 300  20. free  { s20 s21 s69 } } \
                        { n22 550 300  20. free  { s21 s22 s70 } } \
                        { n23 575 300  20. free  { s22 s23 s71 } } \
                        { n24 600 300  20. free  { s23 s24 s72 } } \
                        { n25 625 300  20. fixed { s24     s73 } } \

                        { n26  25 200  20. fixed {     s25 s49 } } \
                        { n27  50 200  20. free  { s25 s26 s50 } } \
                        { n28  75 200  20. free  { s26 s27 s51 } } \
                        { n29 100 200  20. free  { s27 s28 s52 } } \
                        { n30 125 200  20. free  { s28 s29 s53 } } \
                        { n31 150 200  20. free  { s29 s30 s54 } } \
                        { n32 175 200  20. free  { s30 s31 s55 } } \
                        { n33 200 200  20. free  { s31 s32 s56 } } \
                        { n34 225 200  20. free  { s32 s33 s57 } } \
                        { n35 250 200  20. free  { s33 s34 s58 } } \
                        { n36 275 200  20. free  { s34 s35 s59 } } \
                        { n37 300 200  20. free  { s35 s36 s60 } } \
                        { n38 325 200  20. free  { s36 s37 s61 } } \
                        { n39 350 200  20. free  { s37 s38 s62 } } \
                        { n40 375 200  20. free  { s38 s39 s63 } } \
                        { n41 400 200  20. free  { s39 s40 s64 } } \
                        { n42 425 200  20. free  { s40 s41 s65 } } \
                        { n43 450 200  20. free  { s41 s42 s66 } } \
                        { n44 475 200  20. free  { s42 s43 s67 } } \
                        { n45 500 200  20. free  { s43 s44 s68 } } \
                        { n46 525 200  20. free  { s44 s45 s69 } } \
                        { n47 550 200  20. free  { s45 s46 s70 } } \
                        { n48 575 200  20. free  { s46 s47 s71 } } \
                        { n49 600 200  20. free  { s47 s48 s72 } } \
                        { n50 625 200  20. fixed { s48     s73 } } \

                    }

     set spring_def { \
                          { s1  90.0  45.0  10. {n1   n2 } } \
                          { s2  90.0  45.0  10. {n2   n3 } } \
                          { s3  90.0  45.0  10. {n3   n4 } } \
                          { s4  90.0  45.0  10. {n4   n5 } } \
                          { s5  90.0  45.0  10. {n5   n6 } } \
                          { s6  90.0  45.0  10. {n6   n7 } } \
                          { s7  90.0  45.0  10. {n7   n8 } } \
                          { s8  90.0  45.0  10. {n8   n9 } } \
                          { s9  90.0  45.0  10. {n9   n10} } \
                          { s10 90.0  45.0  10. {n10  n11} } \
                          { s11 90.0  45.0  10. {n11  n12} } \
                          { s12 90.0  45.0  10. {n12  n13} } \
                          { s13 90.0  45.0  10. {n13  n14} } \
                          { s14 90.0  45.0  10. {n14  n15} } \
                          { s15 90.0  45.0  10. {n15  n16} } \
                          { s16 90.0  45.0  10. {n16  n17} } \
                          { s17 90.0  45.0  10. {n17  n18} } \
                          { s18 90.0  45.0  10. {n18  n19} } \
                          { s19 90.0  45.0  10. {n19  n20} } \
                          { s20 90.0  45.0  10. {n20  n21} } \
                          { s21 90.0  45.0  10. {n21  n22} } \
                          { s22 90.0  45.0  10. {n22  n23} } \
                          { s23 90.0  45.0  10. {n23  n24} } \
                          { s24 90.0  45.0  10. {n24  n25} } \


                          { s25 90.0  45.0  10. {n26  n27} } \
                          { s26 90.0  45.0  10. {n27  n28} } \
                          { s27 90.0  45.0  10. {n28  n29} } \
                          { s28 90.0  45.0  10. {n29  n30} } \
                          { s29 90.0  45.0  10. {n30  n31} } \
                          { s30 90.0  45.0  10. {n31  n32} } \
                          { s31 90.0  45.0  10. {n32  n33} } \
                          { s32 90.0  45.0  10. {n33  n34} } \
                          { s33 90.0  45.0  10. {n34  n35} } \
                          { s34 90.0  45.0  10. {n35  n36} } \
                          { s35 90.0  45.0  10. {n36  n37} } \
                          { s36 90.0  45.0  10. {n37  n38} } \
                          { s37 90.0  45.0  10. {n38  n39} } \
                          { s38 90.0  45.0  10. {n39  n40} } \
                          { s39 90.0  45.0  10. {n40  n41} } \
                          { s40 90.0  45.0  10. {n41  n42} } \
                          { s41 90.0  45.0  10. {n42  n43} } \
                          { s42 90.0  45.0  10. {n43  n44} } \
                          { s43 90.0  45.0  10. {n44  n45} } \
                          { s44 90.0  45.0  10. {n45  n46} } \
                          { s45 90.0  45.0  10. {n46  n47} } \
                          { s46 90.0  45.0  10. {n47  n48} } \
                          { s47 90.0  45.0  10. {n48  n49} } \
                          { s48 90.0  45.0  10. {n49  n50} } \

                          { s49 90.0  45.0 100. {n1   n26} } \
                          { s50 90.0  45.0  86. {n2   n27} } \
                          { s51 90.0  45.0  84. {n3   n28} } \
                          { s52 90.0  45.0  85. {n4   n29} } \
                          { s53 90.0  45.0 102. {n5   n30} } \
                          { s54 90.0  45.0 115. {n6   n31} } \
                          { s55 90.0  45.0 132. {n7   n32} } \
                          { s56 90.0  45.0 165. {n8   n33} } \
                          { s57 90.0  45.0 190. {n9   n34} } \
                          { s58 90.0  45.0 174. {n10  n35} } \
                          { s59 90.0  45.0 160. {n11  n36} } \
                          { s60 90.0  45.0 153. {n12  n37} } \
                          { s61 90.0  45.0 145. {n13  n38} } \
                          { s62 90.0  45.0 153. {n14  n39} } \
                          { s63 90.0  45.0 160. {n15  n40} } \
                          { s64 90.0  45.0 174. {n16  n41} } \
                          { s65 90.0  45.0 190. {n17  n42} } \
                          { s66 90.0  45.0 165. {n18  n43} } \
                          { s67 90.0  45.0 132. {n19  n44} } \
                          { s68 90.0  45.0 115. {n20  n45} } \
                          { s69 90.0  45.0 102. {n21  n46} } \
                          { s70 90.0  45.0  85. {n22  n47} } \
                          { s71 90.0  45.0  84. {n23  n48} } \
                          { s72 90.0  45.0  86. {n24  n49} } \
                          { s73 90.0  45.0 100. {n25  n50} } \
                      }


     set State(node_list)   {}
     set State(spring_list) {}

     foreach s $spring_def {
         foreach {spring k damp rest nlist} $s {
             set State($spring,name)    $spring
             set State($spring,k)       $k
             set State($spring,kd)      $damp
             set State($spring,rest)    $rest
             set State($spring,nlist)   $nlist
             set State($spring,x1)      0.
             set State($spring,y1)      0.
             set State($spring,x2)      0.
             set State($spring,y2)      0.
             set State($spring,dist)    0.
             set State($spring,force)   0.
             foreach {n1 n2} $nlist {
                 set State($spring,n1)  $n1
                 set State($spring,n2)  $n2
             }
             lappend State(spring_list) $spring
         }
     }

     foreach n $node_def {
         foreach {node x y m type slist} $n {
             set State($node,name)  $node
             set State($node,x)     $x
             set State($node,y)     $y
             set State($node,rad)   [expr $m * 0.1]
             set State($node,mass)  $m
             set State($node,type)  $type
             set State($node,vx)    0.0
             set State($node,vy)    0.0
             set State($node,ax)    0.0
             set State($node,ay)    0.0
             set State($node,Fx)    0.0
             set State($node,Fy)    0.0
             set State($node,slist) $slist
             lappend State(node_list) $node
         }
     }

     foreach s $State(spring_list) {
         set n1 $State($s,n1)
         set n2 $State($s,n2)
         set State($s,x1) $State($n1,x)
         set State($s,y1) $State($n1,y)
         set State($s,x2) $State($n2,x)
         set State($s,y2) $State($n2,y)
     }

 }

 if 0 {

This routine is the main solver. The main loop runs until maxTime is reached. Or until some other exit condition is reached.

At each time step the current forces due to each spring (and gravity) are calculated. A Runge-Kutta routine is then used to calculate the newly calculated positions etc. at the next time step.

 }

 proc solve {func} {
     global State

     if { $State(running) == 0 } {
         return
     }

     set funcq    [uplevel 1 namespace which -command $func]

     set done 0

     while { $done == 0 } {

         set t        $State(t)
         set dmax 0.0

         update_spring_forces
         apply_spring_force_to_nodes

         foreach n $State(node_list) {

             set xvec [list $State($n,vx) $State($n,vy) $State($n,ax) $State($n,ay)]
             set xresult [rungeKuttaStep $n $State(t) $State(dt) $xvec rkfunc]
             foreach {nvx nvy xa ya} $xresult {break}

             set yvec [list $State($n,x) $State($n,y) $nvx $nvy]
             set yresult [rungeKuttaStep $n $State(t) $State(dt) $yvec rkfunc ]
             foreach {nx ny nvx nvy} $yresult {break}

             set m [apply_forces $n $nx $ny $nvx $nvy]

             if {[expr {abs($m)}] > $dmax} {
                 set dmax [expr {abs($m)}]
             }
         }

         set State(dmax) $dmax

         set State(t) [expr {$State(t) + $State(dt)}]

         update_geom

         if {$dmax < $State(tol) || $State(t) > $State(MaxTime)} {
             set done 1
         }
         if { $State(running) == 0 } {
             set done 1
         }
     }
 }

 if 0 {

Calculate the derivatives

 }

 proc rkfunc {n t xvec} {

     foreach {x y vx vy} $xvec {break}
     return [list $vx $vy 0.0 0.0]
 }



 if 0 {

This routine loops over each spring and calculates the force due to the displacement of each node.

 }

 proc update_spring_forces { } {
     global State

     zero_node_forces

     foreach s $State(spring_list) {

         set dx     [expr { $State($s,x2) - $State($s,x1) } ]
         set dy     [expr { $State($s,y2) - $State($s,y1) } ]
         set dist   [expr { sqrt($dx * $dx + $dy * $dy) } ]

         set State($s,dist)  $dist
         set State($s,force) [expr { -$State($s,k) * ($dist - $State($s,rest))}]
     }
 }

 if 0 {

Set the forces applied to each node to zero.

 }

 proc zero_node_forces { } {
     global State

     foreach n $State(node_list) {
         set State($n,Fx) 0.0
         set State($n,Fy) 0.0
     }
 }


 if 0 {

Calculate the forces,accelerations,velocities,etc. to be applied to each node.

 }

 proc apply_spring_force_to_nodes { } {
     global State

     foreach n $State(node_list) {

         set State($n,Fy) [expr {$State(gravity) * $State($n,mass)}]

         foreach s $State($n,slist) {

             if {$n == $State($s,n1)} {
                 set n1 $State($s,n1)
                 set n2 $State($s,n2)
             } else {
                 set n1 $State($s,n2)
                 set n2 $State($s,n1)
             }

             set dx    [expr { $State($n1,x) - $State($n2,x) } ]
             set dy    [expr { $State($n1,y) - $State($n2,y) } ]
             set dvx   [expr { $State($n1,vx) - $State($n2,vx) } ]
             set dvy   [expr { $State($n1,vy) - $State($n2,vy) } ]
             set VdotL [expr { $dvx * $dx + $dvy * $dy } ]

             set dist  $State($s,dist)
             set damp  [expr { $State($s,kd) * $VdotL / $dist }]
             set force [expr { $State($s,force) - $damp } ]

             set Fx [expr { $force * $dx / $dist } ]
             set Fy [expr { $force * $dy / $dist} ]

             set State($n,Fx) [expr {$State($n,Fx) + $Fx}]
             set State($n,Fy) [expr {$State($n,Fy) + $Fy}]
         }
         set State($n,ax) [expr {$State($n,Fx) / $State($n,mass)}]
         set State($n,ay) [expr {$State($n,Fy) / $State($n,mass)}]
     }
 }


 if 0 {

I borrowed this routine from tcllib. Specifically I borrowed the math::calculus::rungeKuttaStep routine so that I could modify it slightly. In the process, I believe that I found a "typo/bug" in the tcllib routine (search MBS below for more info). The bug has been reported and will be fixed in the next release of tcllib (Note: It was in version 1.7).

 }

 ########################################################################

 # rungeKuttaStep --
 #    Integrate a system of ordinary differential equations of the type
 #    x' = f(x,t), where x is a vector of quantities. Integration is
 #    done over a single step according to Runge-Kutta 4th order.
 #
 # Arguments:
 #    t           Start value of independent variable (time for instance)
 #    tstep       Step size of interval
 #    xvec        Vector of dependent values at the start
 #    func        Function taking the arguments t and xvec to return
 #                the derivative of each dependent variable.
 # Return value:
 #    List of values at the end of the step
 #
 proc rungeKuttaStep { n t tstep xvec func } {

    set funcq    [uplevel 1 namespace which -command $func]

    #
    # Four steps:
    # - k1 = tstep*func(t,x0)
    # - k2 = tstep*func(t+0.5*tstep,x0+0.5*k1)
    # - k3 = tstep*func(t+0.5*tstep,x0+0.5*k2)
    # - k4 = tstep*func(t+    tstep,x0+    k3)
    # - x1 = x0 + (k1+2*k2+2*k3+k4)/6
    #
    set tstep2   [expr {$tstep/2.0}]
    set tstep6   [expr {$tstep/6.0}]

    set xk1      [$funcq $n $t $xvec]
    set xvec2    {}
    foreach x1 $xvec xv $xk1 {
       lappend xvec2 [expr {$x1+$tstep2*$xv}]
    }

    set xk2      [$funcq $n [expr {$t+$tstep2}] $xvec2]
    set xvec3    {}
    foreach x1 $xvec xv $xk2 {
       lappend xvec3 [expr {$x1+$tstep2*$xv}]
    }

    set xk3      [$funcq $n [expr {$t+$tstep2}] $xvec3]
    set xvec4    {}
    foreach x1 $xvec xv $xk3 {
       lappend xvec4 [expr {$x1+$tstep *$xv}]
    }
    #***************************************************************************
    # MBS :
    # Previously the above line had:
    #   lappend xvec4 [expr {$x1+$tstep2*$xv}]
    # The "bug" is that $tstep2 in the above line should actually be $tstep
    # From the description above
    #           # - k4 = tstep*func(t+    tstep,x0+    k3)
    # tstep is used rather than tstep*0.5 (ie tstep2)
    #***************************************************************************


    set xk4      [$funcq $n [expr {$t+$tstep}] $xvec4]
    set result   {}
    foreach x0 $xvec k1 $xk1 k2 $xk2 k3 $xk3 k4 $xk4 {
       set dx [expr {$k1+2.0*$k2+2.0*$k3+$k4}]
       lappend result [expr {$x0+$dx*$tstep6}]
    }
    return $result
 }

 if 0 {

Apply the newly calculated positions and velocities of each node.

 }

 proc apply_forces {n nx ny nvx nvy} {
     global State

     set State($n,dx) 0.0
     set State($n,dy) 0.0

     if { $State($n,type) == "fixed" } {
         return 0.0
     }

     set dx [expr { $nx - $State($n,x) } ]
     set dy [expr { $ny - $State($n,y) } ]

     set State($n,dx) $dx
     set State($n,dy) $dy

     set State($n,x) $nx
     set State($n,y) $ny

     set State($n,vx) $nvx
     set State($n,vy) $nvy

     foreach s $State($n,slist) {
         if {$n == $State($s,n1)} {
             set State($s,x1)     $State($n,x)
             set State($s,y1)     $State($n,y)
         } else {
             set State($s,x2)     $State($n,x)
             set State($s,y2)     $State($n,y)
         }
     }
     return [expr {sqrt($dx*$dx + $dy*$dy)} ]
 }

 if 0 {

Initialize the interface :

Hopefully this is pretty straight forward...

 }

 package require Tk
 proc init_display { } {
     global State

     # Window things
     wm title . "TclSpringies"
     wm geometry . "+0+0"
     wm protocol . WM_DELETE_WINDOW {catch "file_quit"}

     # Create canvas
     set State(canvas) [canvas .c -width 700 -height 700 -background black]
     pack $State(canvas) -fill both -expand true

     bind $State(canvas) <ButtonPress-1>   {c_select %W %x %y}
     bind $State(canvas) <B1-Motion>       {c_drag   %W %x %y}
     bind $State(canvas) <ButtonRelease-1> {c_drop   %W %x %y}

     frame .fr -bd 2 -relief raised
     pack .fr -side bottom -fill x
     button .fr.quit  -text "Quit"  -command {exit}
     button .fr.start -text "Start" -command {set State(running) 1; solve rkfunc}
     button .fr.stop  -text "Stop"  -command {set State(running) 0}
     button .fr.reset -text "Reset" -command {reset}
     label .fr.mtlbl -text "Max Time "
     entry .fr.maxt -textvariable State(MaxTime) -width 5

     label .fr.ltol -text "Tolerance "
     entry .fr.etol -textvariable State(tol) -width 8

     pack .fr.quit  -side left
     pack .fr.start -side left
     pack .fr.stop  -side left
     pack .fr.reset -side left
     pack .fr.mtlbl -side left
     pack .fr.maxt  -side left
     pack .fr.ltol  -side left
     pack .fr.etol  -side left


     frame .f1 -bd 2 -relief raised
     label .f1.tlab -text "Time : "
     label .f1.tval -textvariable State(t) -background white -width 20

     pack .f1.tlab -side left
 #    pack .f1.tval -side left -fill x -expand yes
     pack .f1.tval -side left -fill x 

     frame .f2 -bd 2 -relief raised
     label .f2.dlab -text "dMax : "
     label .f2.dval -textvariable State(dmax) -background white -width 20

     pack .f2.dlab -side left
 #    pack .f2.dval -side left -fill x -expand yes
     pack .f2.dval -side left -fill x

     pack .f1 .f2 -side left -fill x -expand yes

     update
 }

 if 0 {

Display the current geometry :

 }

 proc draw_geom { } {
     global State

     set canvasHeight [winfo height $State(canvas)]
     set canvasWidth  [winfo width $State(canvas)]

     foreach s $State(spring_list) {
         set x1 $State($s,x1)
         set y1 [expr $canvasHeight - $State($s,y1)]
         set x2 $State($s,x2)
         set y2 [expr $canvasHeight - $State($s,y2)]
         set State($s,id) [$State(canvas) create line $x1 $y1 $x2 $y2 \
                               -fill white -tag $s]
     }

     foreach n $State(node_list) {
         set r  $State($n,rad)
         set x1 [expr                 $State($n,x) - $r]
         set y1 [expr $canvasHeight - $State($n,y) - $r]
         set x2 [expr                 $State($n,x) + $r]
         set y2 [expr $canvasHeight - $State($n,y) + $r]
         if { $State($n,type) == "fixed" } {
             set State($n,id) [$State(canvas) create oval $x1 $y1 $x2 $y2 \
                                   -outline white -fill red -tag $n ]
         } else {
             set State($n,id) [$State(canvas) create oval $x1 $y1 $x2 $y2 \
                                   -outline white -fill blue -tag $n ]
         }
     }
     update
 }

 if 0 {

Update the geometry after each time step

 }

 proc update_geom { } {
     global State

     set canvasHeight [winfo height $State(canvas)]
     set canvasWidth  [winfo width $State(canvas)]

     foreach s $State(spring_list) {
         set x1 $State($s,x1)
         set y1 [expr $canvasHeight - $State($s,y1)]
         set x2 $State($s,x2)
         set y2 [expr $canvasHeight - $State($s,y2)]
         $State(canvas) coords $State($s,id) $x1 $y1 $x2 $y2 
     }

     foreach n $State(node_list) {
         set dy [expr -1.0 * $State($n,dy)]
         $State(canvas) move $State($n,id) $State($n,dx) $dy
     }
     update
 }

 if 0 {

Some canvas related routines...

Allows the user to select / drag / drop a node using Button-One

 }

 proc c_closest_node {w x y} {
     global State

     set x [$w canvasx $x]
     set y [$w canvasy $y]

     set all_list [$w find all]

     set mdist 9999999.
     set pt_num {}

     foreach item $all_list {

         set ptr  [$w itemcget $item -tags]
         set lndx [lsearch $ptr "n*"]

         if {$lndx != -1} {
             set cord [$w coords $item]
             set nx [expr ([lindex $cord 0] + [lindex $cord 2]) / 2.0]
             set ny [expr ([lindex $cord 1] + [lindex $cord 3]) / 2.0]
             set dx [expr $x - $nx]
             set dy [expr $y - $ny]
             set dst [expr sqrt($dx*$dx + $dy*$dy)]
             if {$dst < $mdist} {
                 set it    $ptr
                 set mdist $dst
             }
         }
     }
     set node [lindex $it 0]
     return $node
 }

 proc c_select {w x y} {
     global State

     set canvasHeight [winfo height $State(canvas)]
     set canvasWidth  [winfo width $State(canvas)]

     set node [c_closest_node $w $x $y]

     set State(curNode) $node
     if {$State(curNode) != "none"} {
         c_move_node $State(curNode) $x $y
     } else {
         set State(curNode) "none"
     }
 }

 proc c_drag {w x y} {
     global State

     set canvasHeight [winfo height $State(canvas)]
     set canvasWidth  [winfo width $State(canvas)]

     if {$State(curNode) != "none"} {
         c_move_node $State(curNode) $x $y
     }
 }

 proc c_drop {w x y} {
     global State

     set canvasHeight [winfo height $State(canvas)]
     set canvasWidth  [winfo width $State(canvas)]

     if {$State(curNode) != "none"} {
         c_move_node $State(curNode) $x $y
     }
     set State(curNode) "none"

 }

 proc c_move_node {n x y} {
     global State

     set canvasHeight [winfo height $State(canvas)]
     set canvasWidth  [winfo width $State(canvas)]

     set r           $State($n,rad)
     set State($n,x) $x
     set State($n,y) [expr {$canvasHeight - $y}]
     set x1 [expr                 $State($n,x) - $r]
     set y1 [expr $canvasHeight - $State($n,y) - $r]
     set x2 [expr                 $State($n,x) + $r]
     set y2 [expr $canvasHeight - $State($n,y) + $r]
     $State(canvas) coords $State($n,id) $x1 $y1 $x2 $y2 

     foreach s $State($n,slist) {
         if {$n == $State($s,n1)} {
             set scoords [$State(canvas) coords $State($s,id)]
             foreach {x1 y1 x2 y2} $scoords {break}
             set State($s,x1)     $x
             set State($s,y1)     [expr {$canvasHeight - $y}]
             $State(canvas) coords $State($s,id) $x $y $x2 $y2
         } else {
             set scoords [$State(canvas) coords $State($s,id)]
             foreach {x1 y1 x2 y2} $scoords {break}
             set State($s,x2)     $x
             set State($s,y2)     [expr {$canvasHeight - $y}]
             $State(canvas) coords $State($s,id) $x1 $y1 $x $y
         }
     }

 }

 if 0 {

Reset to the initial geometry.

 }

 proc reset { } {
     global State

     set State(t) 0.0
     init_geom

     $State(canvas) delete all
     draw_geom
     update_geom
 }

 if 0 {

The 'main' code to start the simulation.

 }

 ########################################################################

 init_display
 init_geom
 draw_geom

 set State(running) 0
 solve rkfunc

 if 0 {

[ Category Command

Category Graphics ]

 }