This is a little toy I wrote for fun. It's a multibody gravity simulation. The most natural use for this is, of course, creating a solar system, so that's what the sample data is. Just a few planets though, no moons so far.
TODO:
using GL or whatever for 3-d would be sweet, but alas, not today.
The function "adddebris" is an attempt to dump small bodies into the L4/L5 points. Doesn't work too well.
In the function "addplanet", distance is in 1e6 km, mass is in 1e24 kg, and period is in earth years.
I hope the real solar system has more accurate numbers or math - the orbits tend to wander a fair amount after a while.
AM You solve the equations using the simplest of all schemes: Euler, first order. If you adapt it Runge-Kutta fourth order (see the Tcllib math::calculus module for an implementation) you will get much better results.
JR I've adapted this to use math::calculus (which interestingly I couldn't find any other uses of on this wiki), so switching beteeen euler and runge-kutta is easy. I don't notice huge differences.
AM: No differences? Surprising. I know that I had, when I did something similar years ago (in FORTRAN using a DOS-based graphics package, so nothing that I would like to reanimate today :). Hm, worth looking into, I think.
As for the use of math::calculus: only few Wiki pages explicitly use external packages. Most simply are standalone. But I welcome the opportunity to see it used here.
package require Tk
package require math::calculus
# gravitational constant
set G .00067
set pi [expr {atan(1)*4}]
proc grav {o1 o2} {
global objs G pi
# objects are list: {X Y mass id dx dy color}
# returns a vector representing the force on o1
foreach {o1x o1y o1m} $o1 {o2x o2y o2m} $o2 {break}
set dx [expr {$o2x-$o1x}]
set dy [expr {$o2y-$o1y}]
set r2 [expr {$dx*$dx+$dy*$dy}]
if ($r2==0) {set r2 1e-10}
set F [expr {($G*($o1m*$o2m)/$r2)/$o1m}]
if {$dx == 0} {
if {$dy > 0} {
set theta [expr {$pi/2}]
} else {
set theta [expr {-$pi/2}]
}
} elseif {$dx < 0} {
set theta [expr {$pi+atan($dy/$dx)}]
} else {
set theta [expr {atan($dy/$dx)}]
}
set Fx [expr {$F*cos($theta)}]
set Fy [expr {$F*sin($theta)}]
return [list $Fx $Fy]
}
proc sumgrav {o} {
# returns the vector sum of all objects on O
global objs
set fv {}
set mo [lindex $objs $o]
foreach obj $objs {
if {[lindex $mo 3] == [lindex $obj 3]} continue
lappend fv [grav $mo $obj]
}
set fx {}
set fy {}
foreach f $fv {
lappend fx [lindex $f 0]
lappend fy [lindex $f 1]
}
return [list [::tcl::mathop::+ {*}$fx] [::tcl::mathop::+ {*}$fy]]
}
proc accelerate {t vec} {
foreach {xv yv xa ya} $vec {break}
return [list $xa $ya 0 0]
}
proc translate {t vec} {
foreach {x y xv yv} $vec {break}
return [list $xv $yv 0 0]
}
set ts 3
set trimdist 5000
set time 0
proc move {} {
global objs ts trimdist time
set oobjs {}
for {set o 0} {$o < [llength $objs]} {incr o} {
set ob [lindex $objs $o]
set tr [sumgrav $o]
# current x, y; mass, name; current x,y velocity;
# x,y acceleration
foreach {cx cy mass name cxv cyv} $ob {xa ya} $tr {break}
set avec [list $cxv $cyv $xa $ya]
set res [math::calculus::rungeKuttaStep $time $ts $avec accelerate]
foreach {nxv nyv xa ya} $res {break}
set tvec [list $cx $cy $nxv $nyv]
set res [math::calculus::rungeKuttaStep $time $ts $tvec translate]
foreach {nx ny nxv nyv} $res {break}
if {sqrt($nx*$nx+$ny*$ny) < $trimdist} {
set ob [concat $nx $ny [lrange $ob 2 3] $nxv $nyv [lrange $ob 6 end]]
lappend oobjs $ob
}
}
set objs $oobjs
}
set scale 0.2
set trails 0
set okeep 1000
proc animate {} {
global trails okeep
if {!$trails} {
.c delete all
} else {
set all [.c find all]
if {[llength $all] > $okeep} {
.c delete {*}[lrange $all 0 [expr {[llength $all]-1000}]]
}
}
global objs scale
set ccx [expr {[winfo width .c]/2}]
set ccy [expr {[winfo height .c]/2}]
foreach o $objs {
set cx [expr {[lindex $o 0]*$scale+$ccx}]
set cy [expr {[lindex $o 1]*$scale+$ccy}]
set cr [expr {[lindex $o 2]*$scale/100000}]
.c create oval [expr {$cx-$cr}] [expr {$cy-$cr}] [expr {$cx+$cr}] [
expr {$cy+$cr}] -outline [lindex $o 6]
}
}
set playing 1
set delay 25
proc playing {} {
global delay
move
animate
global playing
if {$playing} {after $delay [namespace code [info level 0]]}
}
proc float number {
expr {$number+0.0}
}
proc addplanet {name distance mass period color} {
# distance is in 1e6 km
# mass is in 1e24 kg
# period is in earth years
global objs
set speed [expr {$distance*1000*2*3.14/($period*365*846)}]
lappend objs [list 0 [float $distance] [float $mass] $name $speed 0 $color]
}
proc adddebris tag {
global objs
set ob [lindex $objs [lsearch $objs *$tag*]]
foreach {ox oy} $ob {sx sy} [lindex $objs 0] {break}
set d 10
set dx [expr {$ox-$sx}]
set dy [expr {$oy-$sy}]
set th [expr {atan($dy/($dx+1e-10))}]
set dist [expr {sqrt($dx*$dx+$dy*$dy)}]
for {set c 0} {$c < 5} {incr c} {
lappend objs [list [expr {cos($th+1.047)*$dist+$d*rand()}] [
expr {sin($th+1.047)*$dist+$d*rand()}] 1e-10 trojan 0 0 gray]
lappend objs [list [expr {cos($th-1.047)*$dist+$d*rand()}] [
expr {sin($th-1.047)*$dist+$d*rand()}] 1e-10 trojan 0 0 gray]
}
}
proc sysmom {} {
global objs
set mx 0
set my 0
foreach o $objs {
set mx [expr {$mx + [lindex $o 2]*[lindex $o 4]}]
set my [expr {$my + [lindex $o 2]*[lindex $o 5]}]
}
list $mx $my
}
proc recenter {} {
global objs
set sun [lindex $objs 0]
set sx [lindex $sun 0]
set sy [lindex $sun 1]
set nobjs {}
foreach o $objs {
lappend nobjs [concat [expr {[lindex $o 0]-$sx}] [
expr {[lindex $o 1]-$sy}] [lrange $o 2 end]]
}
set objs $nobjs
}
proc setup {} {
canvas .c -background black -width 400 -height 400
frame .panel
pack .c -expand t -fill both
pack .panel
# controls
label .zooml -text "Zoom"
scale .zoom -resolution .005 -orient horizontal \
-from .001 -to 2 -variable scale
grid .zooml .zoom -in .panel -sticky s
label .speedl -text "Delay"
scale .speed -resolution 1 -orient horizontal -from 1 -to 100 -variable delay
grid .speedl .speed -in .panel -sticky s
label .timel -text "Time step"
scale .time -resolution 0.5 -orient horizontal -from 0.5 -to 50 -variable ts
grid .timel .time -in .panel -sticky s
label .playl -text "Animate"
checkbutton .play -command playing -variable playing
grid .playl .play -in .panel
label .traill -text "Show Trails"
checkbutton .trail -variable trails
grid .traill .trail -in .panel
button .clear -command ".c delete all" -text Clear
button .center -command "recenter" -text Recenter
grid .clear .center -in .panel
}
# solar system
# the strange number for the sun's velocity is to try to keep the overall
# momentum of the system to zero
set objs { {0 0 1989000 sun -0.00163355465017 0 yellow} }
addplanet mercury 57.9 .33 .23 red
addplanet venus 108.2 4.86 .613 green
addplanet earth 149 5.97 1 blue
addplanet mars 227 .64 2 red
addplanet jupiter 778 1900.0 11.868 green
addplanet saturn 1429 568.0 29.47 yellow
addplanet uranus 2870 86.8 84.06 blue
addplanet neptune 4504 102.47 164.90 green
addplanet pluto 5913 .012 248.08 grey
setup
playing
schlenk: Hmm, may i suggest OpenGL support for 3D?
TV: Nice, though a bit depressing after a while.. It could be a good idea to have a course screen resolution as a test, I had that when I was a teenager doing this sort of thing, first order Euler or not (or backward, with prediction) then it is clear where roundoffs come from, or (energy) compensate by.
Considering the attraction falls of with a square as a function of distance, it can be considered that using complicated or higher order approximations (unless they are for this purpose) may leave you with more stable (not usually the case with higher order polynomials, though 4 is not that high) or more accurate solutions, they also may not take non-linearity better into account than a lower order solution with a smaller timestep, and in extreme cases, like when planets sheer just past eachother, may be wrong.
I didn't look at the code, but would it be doable to make longer, and maybe fading trails ?
JR: longer trails is a cinch, just change the value of "okeep". Increasing it too much will eventually bog down the canvas widget tho. Fading trails would be somewhat more work and significantly slower (it would need to run through the whole list of objects on the canvas with each or every few passes)
GWM: If you need more accuracy you could try the Runge-Kutta method.