from visual import *

scene.range=5
scene.autoscale=0

star=[]
trail=[]
star.append(sphere(color=color.red,radius=0.2,pos=(-0.5,0,0),m=1))
trail.append(curve(color=color.red,pos=[star[0].pos]))
star[0].p=star[0].m*vector(0,0,0)
star.append(sphere(color=color.blue,radius=0.2,pos=(0.5,0,0),m=0.0001))
trail.append(curve(color=color.blue,pos=[star[1].pos]))
star[1].p=star[1].m*vector(0,1,0)

nstars=2

t=0
dt=0.001
tmax=20

while (t < tmax):
    rate(1./dt)
    for i in range(nstars):
        star[i].F=vector(0,0,0)
        for j in range(nstars):
            if (i != j):
                star[i].F=star[i].F+star[i].m*star[j].m*(star[j].pos-star[i].pos)/pow(mag(star[j].pos-star[i].pos),3)
    for i in range(nstars):
        star[i].pos=star[i].pos+star[i].p/star[i].m*dt
        trail[i].append(star[i].pos)
        star[i].p=star[i].p+star[i].F*dt

                

    t=t+dt