A naive simulator of gravity, in less than 200 Python lines...

#!/usr/bin/env python
"""
A naive simulator of gravity, in less than 200 lines

Coded by Thanassis Tsiodras, during a 2010 summer night...
(http://users.softlab.ntua.gr/~ttsiod)

Well, 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 gravity law, and see
what happens in... a universe of my own making...

Well, here's what I did: 
The following code 'sprays' 200 'particles' randomly, inside a
900x600 window (the values are below, change them at will).
Afterwards, it applies a very simple set of laws:

- Gravity, inversely proportional to the square of the distance,
  and linearly proportional to the product of the two masses
- Elastic collissions of two objects if they are close enough to
  touch: a merged object is then created, that maintains the 
  momentum (mass*velocity) and the mass of the two merged ones.

Well, I probably goofed at something, but I like the results -
they appear very natural to me, and more importantly, after 
half a minute or so, the remaining objects have such momentum,
that they leave the original 'window' with great speeds, travelling
into space...

Yes, I've proven the Big Bang!
I want my Nobel prize now, thank you :-)

(When you reach this stage, use the 'A' and 'Z' keys to zoom
in/out and see the 'expansion of the universe'...)

Released under the GPL, code available at:
http://users.softlab.ntua.gr/~ttsiod/gravity.html

"""

import sys
import math
import pygame
import random

# The window size
WIDTH = 900
HEIGHT = 600
WIDTHD2 = WIDTH/2.
HEIGHTD2 = HEIGHT/2.

# The number of simulated particles
PARTICLES = 200

# The density of the particles - used to calculate their mass 
# from their volume (i.e. via their radius)
DENSITY = 0.001

# The gravity coefficient. 
# It's my universe, I can pick whatever I want.
GRAVITYSTRENGTH = 1.e4

class Particle:
    def __init__(self):
	# Position
	self._x  = float(random.randint(0,WIDTH))
	self._y  = float(random.randint(0,HEIGHT))
	# Speed
	self._vx = 0.  
	self._vy = 0.
	# Radius
	self._radius = 1.5
	# Use density to calculate mass from radius
	self.setMassFromRadius()
	# When two particles are close enough to touch:
	# - _merged will be set to True for the lighter one
	# - it will no longer be simulated, since its mass
	#    will be added to the larger one, and its momentum 
	#    (mass*velocity) will be merged to the larger one's
	#    (see code for merging particles below, in main loop)
	self._merged = False
    def resetAccel(self):
	# Acceleration, set to 0 in each "calculate total force" loop
	self._ax = 0.0
	self._ay = 0.0
    def addAccel(self, particle):
	# Calculate the acceleration contribution from another particle
	dx = particle._x - self._x
	dy = particle._y - self._y
	dsq = dx*dx + dy*dy  # distance squared
	dr = math.sqrt(dsq)  # distance
	force = GRAVITYSTRENGTH*self._mass*particle._mass/dsq if dsq>1e-10 else 0.
	# Accumulate acceleration...
	self._ax += force*dx/dr
	self._ay += force*dy/dr
    def updatePos(self):
	# Afrer "addAccel" has been called to account for all other particles,
	# update my position and my speed...
	self._x += self._vx
	self._y += self._vy
	self._vx += self._ax
	self._vy += self._ay
	self.resetAccel() # prepare for next round
    def setMassFromRadius(self):
	# The volume is (4/3)*Pi*(r^3), so...
	self._mass = DENSITY * 4. * math.pi * (self._radius**3.)/3.
    def setRadiusFromMass(self):
	# Reversing the previous formula, to calculate radius from mass
	# (used after merging of two particles - mass is added, and new
	# radius is calculated from this formula)
	self._radius = (3.*self._mass/(DENSITY * 4.*math.pi))**(0.3333)

def panic(msg):
    sys.stderr.write(msg+"\n")
    sys.exit(1)

def main():
    pygame.init() 
    win=pygame.display.set_mode((WIDTH,HEIGHT))

    # And God said: Let there be lights in the firmament of the heavens...
    particles = []
    keysPressed = {pygame.K_a:0, pygame.K_z:0, pygame.K_ESCAPE:0}
    for i in xrange(0,PARTICLES):
	particles.append(Particle())

    # Zoom factor, changed at runtime via the 'A' and 'Z' keys
    zoom = 1.0

    while True:
	pygame.display.flip()
	win.fill((0,0,0))
	win.lock()
	for particle in particles:
	    if not particle._merged: # for those that have not been merged,
		# Draw a circle as big as their radius, but take zoom 
		# factor into account
		pygame.draw.circle(win, (255,255,255), 
		    (int(WIDTHD2 + zoom*WIDTHD2*(particle._x-WIDTHD2)/WIDTHD2), 
		     int(HEIGHTD2 + zoom*HEIGHTD2*(particle._y-HEIGHTD2)/HEIGHTD2)), 
		     int(particle._radius*zoom), 0)
	win.unlock()
	while True:

	    # Update the keysPressed dictionary, so we know what keys are pressed
	    evt = pygame.event.poll()
	    if evt.type == pygame.NOEVENT:
		break
	    elif evt.type == pygame.KEYDOWN:
		keysPressed[evt.key] = 1
	    elif evt.type == pygame.KEYUP:
		if keysPressed.has_key(evt.key) and keysPressed[evt.key] == 1:
		    keysPressed[evt.key] = 0
		break

	# Update all particle positions and speeds, by running a loop on each one...
	for i in xrange(0, len(particles)):
	    if particles[i]._merged: continue
	    # ...and calculating the contributions of all the others to its velocity
	    # (via the gravity force)
	    particles[i].resetAccel()
	    for j in xrange(0, len(particles)):
		if i == j or particles[j]._merged: continue # not pulled by myself :-)
		particles[i].addAccel(particles[j])
	    particles[i].updatePos()

	# See if we should merge the ones that are close enough to touch...
	for i in xrange(0, len(particles)):
	    if particles[i]._merged: continue
	    for j in xrange(0, len(particles)):
		if i == j or particles[j]._merged: continue
		dx = particles[i]._x - particles[j]._x
		dy = particles[i]._y - particles[j]._y
		dsq = dx*dx + dy*dy
		dr = math.sqrt(dsq)
		if dr<=(particles[i]._radius + particles[j]._radius):
		    # Merge us, we are close enough to touch!
		    if particles[i]._mass > particles[j]._mass:
			ii,jj = i,j
		    else:
			ii,jj = j,i
		    # ii is the index of the biggest one (mass-wise)
		    particles[jj]._merged = True
		    # maintain the momentum - elastic collisions...
		    newvx = ( particles[ii]._vx*particles[ii]._mass + \
			      particles[jj]._vx*particles[jj]._mass) / \
			        (particles[ii]._mass + particles[jj]._mass)
		    newvy = ( particles[ii]._vy*particles[ii]._mass + \
			      particles[jj]._vy*particles[jj]._mass) / \
			        (particles[ii]._mass + particles[jj]._mass)
		    # maintain the mass (just add them)
		    particles[ii]._mass += particles[jj]._mass
		    # based on the new mass, figure out the new radius
		    particles[ii].setRadiusFromMass()
		    # update velocity of merged object, to maintain momentum
		    particles[ii]._vx = newvx
		    particles[ii]._vy = newvy

	# if 'A' or 'Z' are pressed, update zoom factor:
	if keysPressed[pygame.K_a]:
	    zoom /= 0.99
	    print "Zoom factor:", zoom
	if keysPressed[pygame.K_z]:
	    zoom /= 1.01
	    print "Zoom factor:", zoom
	if keysPressed[pygame.K_ESCAPE]:
	    break
	if evt.type == pygame.NOEVENT:
	    pygame.time.wait(20) # don't get beyond 50 frames per second...
	    continue

if __name__ == "__main__":
    try:
        import psyco
        psyco.profile()
    except:
        print 'Psyco not found, ignoring it'
    main()