Gravity

Over on the flickr nodebox group there is a little discussion about a particle library for nodebox. I have played with some very basic gravity simulations and Tom asked if I had any useful code. The answer is: 'probably not,' but I'm going to post some anyway--who knows maybe someone will refine it and do something interesting. This is an animated gravity simulation. It should be self-explanatory and it's very simple.

speed(1550)
size(800,800)
strokewidth(.05)
from Numeric import *
class System:
    '''
    System is a class that keeps track of the forces between Particle Objects
    '''
    def __init__(self):
        self.universe = []
        self.center = None
        self.G = 10
    def addBody(self, *bodies):
        for body in bodies:
            self.universe.append(body)
    def distance(self, body1, body2):
        dist = sqrt(sum((body1.position - body2.position)**2))
        # We pretend collisions don't happen
        # but without collisions distance get artificially small and gravity becomes too large
        if dist < 10: return 10
        else: return dist
    def gravity(self, body1, body2, dist):
        ''' Just like Newton taught us '''
        return (body1.mass * body2.mass) / (dist**2) * self.G
    def makeGravity(self):
        ''' cycle though each pair of bodies in system 
        and determine force--is there a better way?'''
        for i, body1 in enumerate(self.universe):
            for body2 in self.universe[i+1:]:
                dist = self.distance(body1, body2)
                gForce = self.gravity(body1, body2, dist )
                body1.receiveForceVector (gForce * (body1.position - body2.position) /dist)
                body2.receiveForceVector (gForce * (body2.position - body1.position) /dist)
        map(lambda b: b.updatePosition(1), self.universe)
          
class Particle:
    def __init__(self, mass, position, velocity):
        '''velocity and position should be x,y,z vectors. It would probably be a good idea to normalize mass'''
        self.mass = mass
        self.velocity = array(velocity, Float)
        self.position = array(position, Float)
        self.x = self.position[0]
        self.y = self.position[1]
        self.z = self.position[2]
        self.isStationary = False
    def draw(self):
        oval(self.position[0]-2, self.position[1]-2, 4, 4 )
    def receiveForceVector(self, force):
        """ force should be individual vector components i.e [.2 .3 0] """
        self.velocity -= force/self.mass
    def updatePosition(self, t):
        if not self.isStationary:
            self.position += self.velocity * t
            self.x, self.y, self.z = self.position
        self.draw()
       
def setup():        
    global p1_mass, p2_mass, p3_mass, a, p1, p2, p3, p4, system
    p1_mass=.4
    p2_mass=.1
    p3_mass=.1
    a = .45 
    system = System()
    
    system.addBody( Particle(p1_mass, (-300, 0, 0), (0,-a,0)),
                    Particle(p2_mass, (-300, 10, 0), (0, -a,0)), 
                    Particle(p3_mass, (-300, -9, 0), (1.01*a,0, 0))
                    )
    
    sun  = Particle(10, (0, -0, 0), (0, 0, 0))
    sun.isStationary = True
    system.addBody(sun)
 
def draw():
    global p1_mass, p2_mass, p3_mass, a, p1, p2, p3, system
    translate(400, 400)
    nofill()
    stroke(.2, .2, .5, .5)
    system.makeGravity()

Hi Mark,

Thanks for the code! I'm going to have a look at it.
I've also come across PyODE (http://pyode.sourceforge.net/) which might be interesting. Would be nice to have a pure Python implementation though... I still hope to get an OK from Traer to port his physics in Processing so we don't have to reinvent the wheel.