A naive simulator of gravity, written in Python

(Coded during a 2010 summer night... updated with Runge-Kutta solver during Xmas 2012 weekend... Reddit-ed

I've always been fascinated by space - ever since I read "The Family of the Sun", when I was young. And I always wanted to simulate what I've read about Newton's law of gravity, and see what happens in... a universe of my own making... :‑)

First with 30 planets, then with just one (to show stable orbit, thanks to RK4

Gravity, inversely proportional to the square of the distance, and linearly proportional to the product of the two masses

Elastic "mergings" of two planets if they are close enough to touch: a merged planet is then created, that maintains the momentum (mass times velocity) and the mass of the two merged ones.

self._x += self._vx self._y += self._vy self._vx += self._ax self._vy += self._ay

So eventually, I sat down and did it. The following code "sprays" planets randomly, inside a 900x600 window (the values are in the code, change them at will). It then proceeds to apply the following:My original version used Euler integration - i.e. assuming a unit time step......it just applied a crude form of Euler's method in updating position and velocity. I always meant to update the code with a more stable solution... and I finally found some time to hack it, during the Xmas 2012 weekend: the code now uses the RK4 solution of the velocity / acceleration differential equation. I based the patch on the excellent work of Glenn Fiedler (definitely worth checking out, if you like this sort of thing).

And now... I can observe stable orbits - and have long-term fun watching chaotic solar systems calm down into peaceful ones like our own :‑)

Download / usage

gravityRK4.exe 20

sudo apt-get install python-pygame

bash$ ./gravityRK4.py 20

If you use Windows, download a pre-compiled version from here , and run with(20 is the number of planets - change it to as many as you want). Otherwise, the Python code (less than 300 commented lines) is here . After downloading, make sure you have pygame installed (or whatever your distro/environment needs) and then, just......to run a simulation with 20 planets. Change this number up to whatever number your machine can handle (think of it as a benchmark :‑)

You can also use the numeric keypad's +/- to zoom in/out, and press SPACE to toggle showing/hiding the orbits trace.

Enjoy!

(More of my "semi-scientific" models of natural processes can be found here.)





import sys import math import pygame import random from collections import defaultdict WIDTH , HEIGHT = 900 , 600 WIDTHD2 , HEIGHTD2 = WIDTH / 2 ., HEIGHT / 2 . PLANETS = 30 DENSITY = 0.001 GRAVITYSTRENGTH = 1 . e4 g_listOfPlanets = [] class State : def __init__ ( self , x , y , vx , vy ): self . _x , self . _y , self . _vx , self . _vy = x , y , vx , vy def __repr__ ( self ): return 'x:{x} y:{y} vx:{vx} vy:{vy}' . format ( x = self . _x , y = self . _y , vx = self . _vx , vy = self . _vy ) class Derivative : def __init__ ( self , dx , dy , dvx , dvy ): self . _dx , self . _dy , self . _dvx , self . _dvy = dx , dy , dvx , dvy def __repr__ ( self ): return 'dx:{dx} dy:{dy} dvx:{dvx} dvy:{dvy}' . format ( dx = self . _dx , dy = self . _dy , dvx = self . _dvx , dvy = self . _dvy ) class Planet : def __init__ ( self ): if PLANETS == 1 : self . _st = State ( 150 , 300 , 0 , 2 ) else : self . _st = State ( float ( random . randint ( 0 , WIDTH )), float ( random . randint ( 0 , HEIGHT )), float ( random . randint ( 0 , 300 )/ 100 .)- 1.5 , float ( random . randint ( 0 , 300 )/ 100 .)- 1.5 ) self . _r = 1.5 self . setMassFromRadius () self . _merged = False def __repr__ ( self ): return repr ( self . _st ) def acceleration ( self , state , unused_t ): ax = 0.0 ay = 0.0 for p in g_listOfPlanets : if p is self or p . _merged : continue dx = p . _st . _x - state . _x dy = p . _st . _y - state . _y dsq = dx * dx + dy * dy dr = math . sqrt ( dsq ) force = GRAVITYSTRENGTH * self . _m * p . _m / dsq if dsq > 1e-10 else 0 . ax += force * dx / dr ay += force * dy / dr return ( ax , ay ) def initialDerivative ( self , state , t ): ax , ay = self . acceleration ( state , t ) return Derivative ( state . _vx , state . _vy , ax , ay ) def nextDerivative ( self , initialState , derivative , t , dt ): state = State ( 0 ., 0 ., 0 ., 0 .) state . _x = initialState . _x + derivative . _dx * dt state . _y = initialState . _y + derivative . _dy * dt state . _vx = initialState . _vx + derivative . _dvx * dt state . _vy = initialState . _vy + derivative . _dvy * dt ax , ay = self . acceleration ( state , t + dt ) return Derivative ( state . _vx , state . _vy , ax , ay ) def updatePlanet ( self , t , dt ): a = self . initialDerivative ( self . _st , t ) b = self . nextDerivative ( self . _st , a , t , dt * 0.5 ) c = self . nextDerivative ( self . _st , b , t , dt * 0.5 ) d = self . nextDerivative ( self . _st , c , t , dt ) dxdt = 1.0 / 6.0 * ( a . _dx + 2.0 *( b . _dx + c . _dx ) + d . _dx ) dydt = 1.0 / 6.0 * ( a . _dy + 2.0 *( b . _dy + c . _dy ) + d . _dy ) dvxdt = 1.0 / 6.0 * ( a . _dvx + 2.0 *( b . _dvx + c . _dvx ) + d . _dvx ) dvydt = 1.0 / 6.0 * ( a . _dvy + 2.0 *( b . _dvy + c . _dvy ) + d . _dvy ) self . _st . _x += dxdt * dt self . _st . _y += dydt * dt self . _st . _vx += dvxdt * dt self . _st . _vy += dvydt * dt def setMassFromRadius ( self ): self . _m = DENSITY * 4 .* math . pi *( self . _r ** 3 .)/ 3 . def setRadiusFromMass ( self ): self . _r = ( 3 .* self . _m /( DENSITY * 4 .* math . pi ))**( 0.3333 ) def main (): pygame . init () win = pygame . display . set_mode (( WIDTH , HEIGHT )) keysPressed = defaultdict ( bool ) def ScanKeyboard (): while True : evt = pygame . event . poll () if evt . type == pygame . NOEVENT : break elif evt . type in [ pygame . KEYDOWN , pygame . KEYUP ]: keysPressed [ evt . key ] = evt . type == pygame . KEYDOWN global g_listOfPlanets , PLANETS if len ( sys . argv ) == 2 : PLANETS = int ( sys . argv [ 1 ]) g_listOfPlanets = [] for i in xrange ( 0 , PLANETS ): g_listOfPlanets . append ( Planet ()) def planetsTouch ( p1 , p2 ): dx = p1 . _st . _x - p2 . _st . _x dy = p1 . _st . _y - p2 . _st . _y dsq = dx * dx + dy * dy dr = math . sqrt ( dsq ) return dr <=( p1 . _r + p2 . _r ) sun = Planet () sun . _st . _x , sun . _st . _y = WIDTHD2 , HEIGHTD2 sun . _st . _vx = sun . _st . _vy = 0 . sun . _m *= 1000 sun . setRadiusFromMass () g_listOfPlanets . append ( sun ) for p in g_listOfPlanets : if p is sun : continue if planetsTouch ( p , sun ): p . _merged = True zoom = 1.0 t , dt = 0 ., 1 . bClearScreen = True pygame . display . set_caption ( 'Gravity simulation (SPACE: show orbits, ' 'keypad +/- : zoom in/out)' ) while True : t += dt pygame . display . flip () if bClearScreen : win . fill (( 0 , 0 , 0 )) win . lock () for p in g_listOfPlanets : if not p . _merged : pygame . draw . circle ( win , ( 255 , 255 , 255 ), ( int ( WIDTHD2 + zoom * WIDTHD2 *( p . _st . _x - WIDTHD2 )/ WIDTHD2 ), int ( HEIGHTD2 + zoom * HEIGHTD2 *( p . _st . _y - HEIGHTD2 )/ HEIGHTD2 )), int ( p . _r * zoom ), 0 ) win . unlock () ScanKeyboard () for p in g_listOfPlanets : if p . _merged or p is sun : continue p . updatePlanet ( t , dt ) for p1 in g_listOfPlanets : if p1 . _merged : continue for p2 in g_listOfPlanets : if p1 is p2 or p2 . _merged : continue if planetsTouch ( p1 , p2 ): if p1 . _m < p2 . _m : p1 , p2 = p2 , p1 p2 . _merged = True if p1 is sun : continue newvx = ( p1 . _st . _vx * p1 . _m + p2 . _st . _vx * p2 . _m )/( p1 . _m + p2 . _m ) newvy = ( p1 . _st . _vy * p1 . _m + p2 . _st . _vy * p2 . _m )/( p1 . _m + p2 . _m ) p1 . _m += p2 . _m p1 . setRadiusFromMass () p1 . _st . _vx , p1 . _st . _vy = newvx , newvy if keysPressed [ pygame . K_KP_PLUS ]: zoom /= 0.99 if keysPressed [ pygame . K_KP_MINUS ]: zoom /= 1.01 if keysPressed [ pygame . K_ESCAPE ]: break if keysPressed [ pygame . K_SPACE ]: while keysPressed [ pygame . K_SPACE ]: ScanKeyboard () bClearScreen = not bClearScreen verb = "show" if bClearScreen else "hide" pygame . display . set_caption ( '%s orbits, keypad +/- : zoom in/out)' % verb ) if __name__ == "__main__" : try : import psyco psyco . profile () except : print 'Psyco not found, ignoring it' main ()

Back to index My CV Last update on: Sat Jul 20 15:34:11 2019

The comments on this website require the use of JavaScript. Perhaps your browser isn't JavaScript capable or the script is not being run for another reason. If you're interested in reading the comments or leaving a comment behind please try again with a different browser or from a different connection.