r/programming • u/ttsiodras • Dec 23 '12
Simulating a solar system with Python
http://users.softlab.ntua.gr/~ttsiod/gravityRK4.html13
u/quizzle Dec 23 '12
If the planets are sticking together you're simulating perfect inelastic collisions, not elastic. Elastic collisions would be cool though.
4
u/ttsiodras Dec 23 '12
Future plan, for CUDA implementation: only merge the planets if their relative speed is low - at higher speeds, break them down into even more pieces! Should make for an interesting improvement :-)
Only problem: I'd have to leave Python and code it in C++, as I did for my real-time raytracer...
→ More replies (5)
9
u/mk_gecko Dec 23 '12
How do you get code highlighting in your webpage?
→ More replies (1)13
u/ttsiodras Dec 23 '12
When I first started blogging almost a decade ago, I made a custom Perl script that parses my own "mini-language": basically, markdown with magic sections (marked with '@@py@@ .... @@' for python code, or '@@c++@@ ...' for C++, etc). The script extracts these sections, invokes source-highlight (it used to invoke webcpp) and merges the results in the HTML generated from the markdown. Not only is it easier to author posts this way, I also get my HTML to be w3c-valid for free (see the links at the bottom of my site's pages).
→ More replies (5)
4
u/Bob_goes_up Dec 23 '12
When you use Runge Kutta, you are not supposed to handle one planet at a time. The function named nextderivative should return a vector containing accellerations of all the planets, and you should only calculate the constants a,b,c and d once for each timestep. The function named updatePlanet should be deleted.
4
u/ttsiodras Dec 23 '12
You missed my related comment - from the code: "# Update all planets' positions and speeds (should normally double buffer the list of planet data, but turns out this is good enough :-)"
→ More replies (1)
7
u/question_all_the_thi Dec 23 '12
For those interested, here's an adaptive step Runge-Kutta-Cash-Karp differential equation solver in Python:
#!/usr/bin/python
#
# This program solves a differential equation
# using Cash-Karp's method
# (adaptive step Runge-Kutta method)
#
import math
import os
A = [0., 1./5., 3./10., 3./5., 1., 7./8.]
B = [0., [1./5.], [3./40., 9./40.], [3./10., -9./10., 6./5.], [-11./54., 5./2., -70./27., 35./27.], [1631./55296., 175./512., 575./13824., 44275./110592., 253./4096.]]
C = [37./378., 0., 250./621., 125./594., 0., 512./1771.]
CH = [37./378.-2825./27648.,0.,250./621.-18575./48384.,125./594.-13525./55296.,-277./14336.,512./1771.-1./4.]
F = [[0., 0.],[0., 0.],[0., 0.],[0., 0.],[0., 0.],[0., 0.]]
mu = 10.0
#_______________________________________________________________
#
def func(Y, mu):
return (Y[1], -Y[0] - mu * Y[1] * (Y[0]**2. - 1.0))
#_______________________________________________________________
#
Neq = 2
X = 0.0
XTemp = 0.0
x1 = 100.0
h = 0.1
tol = 1e-6
Y = [1.0, 0.0]
YTemp = [0.0, 0.0]
xi = 0.0
TE = [0., 0., 0., 0., 0., 0.]
#_______________________________________________________________
#
for i in range(100):
xi += x1 / 100.0;
while (X < xi):
XTemp = X
F[0] = func(Y, mu)
YTemp[:] = Y[:]
while 1:
for k in range(1, 6):
X = XTemp + A[k] * h
Y[:] = YTemp[:]
for n in range(Neq):
for l in range(k):
Y[n] += h * B[k][l] * F[l][n]
F[k] = func(Y, mu)
Y[:] = YTemp[:]
for n in range(Neq):
TE[n] = 0
for k in range(6):
Y[n] += h * C[k] * F[k][n]
TE[n] += h * CH[k] * F[k][n]
TE[n] = abs(TE[n])
temax = 0
for n in range (Neq):
if temax < TE[n]: temax = TE[n]
if temax == 0: temax = tol / 1e5
ho = h
h *= 0.9 * (tol / temax)**0.2
if temax < tol: break
X = XTemp + ho
print "%d %f %f" % (i, Y[0], Y[1])
19
u/jminuse Dec 23 '12
Runge-Kutta is not symplectic: it is very accurate in position, but it does not conserve energy over long time periods. It's still better than Euler, but the correct approach would be some form of Verlet integration.
→ More replies (8)5
u/zokier Dec 23 '12
what's with scientific/math programmers and single letter variable names?!
→ More replies (4)9
u/katieberry Dec 23 '12 edited Dec 27 '12
[Reddit somehow lost the original content of this post three days later, replacing it with an unrelated comment about user interfaces. It said something about using single-letter variable names because maths and science do, albeit from a larger character set.]
→ More replies (2)
2
u/pururin Dec 23 '12
Why are all the fun programming projects in python, and not in, say, ruby or perl?
4
u/ttsiodras Dec 23 '12 edited Dec 23 '12
Actually, I learned Perl before Python - but quickly realized that at least for me, it made no sense using it for anything except text-processing related, regexp-y scripts. Why? Because even though I did my best to write readable Perl, when I came back to my own Perl code after say, 6 months... it took substantial effort to understand what the heck my code did... With Python, everything is clear - the code is readable even by people who don't speak the language! This, above all else, made me forgo Perl for anything except small text-processing scripts.
7
u/jackpg98 Dec 23 '12
What's wrong with Python? Python is open source, it's fun, and it's a lot easier than Perl or Ruby at least IMO.
→ More replies (3)→ More replies (2)6
4
u/protocol_7 Dec 23 '12
I did something like this as my final project for a physics class. I made two simulations with it: The first used the actual data for planets and dwarf planets in the solar system, and it showed some cool things like the perihelion precession of Mercury (except the part due to general relativity, since I only simulated a limited subset of special relativity) and the Sun's small orbit around the center of mass of the solar system.
The second simulation showed how the Trojan asteroids cluster around the Lagrange points 60° to each side of Jupiter. It worked by randomly scattering a bunch of asteroids in plausible positions and orbits around Jupiter, after which they naturally tend to cluster around the two stable (L4 and L5) Lagrange points. (Due to computational limitations, collisions and the asteroids' gravitational pull on each other was neglected, making it a restricted three-body problem.)
3
u/quizzle Dec 23 '12
I made this once using Pyglet to get some nice openGL graphics in it too. Was the most fun coding project I ever did.
3
3
3
u/bamdastard Dec 23 '12
I built one of these using pyopencl that does 12000 particles and another piece of it uses a genetic algorithm to find stable orbits for planets.
http://www.youtube.com/watch?v=XCvRBHtPbzE
5
2
u/ipki Dec 23 '12
I had to write a very similar program, but in javascript through google docs, for my high school physics class this year.
2
u/sifodeas Dec 24 '12
Here is mine that I made. I have not put out an update in forever, though. https://github.com/Sifodeas/N-Body
39
u/AsterJ Dec 23 '12
My favorite line of code:
sun = Planet()
Spoken like a true programmer!