GlowScript 2.7 VPython

scene.range=2e12

G = 6.7e-11 

bigstar = sphere(pos=vector(-5e11,0,0), radius=4e10, color=color.red,make_trail=False, interval=10)
bigstar.mass = 3e30
bigstar.p = vector(0, 3e3, 0) * bigstar.mass

smallstar = sphere(pos=vector(5e11,0,0), radius=3e10, color=color.yellow, make_trail=False, interval=10)
smallstar.mass = 1e30
smallstar.p = vector(0, -9e3, 0) * smallstar.mass

asteroid = sphere(pos=vector(1.43e12,0,0), radius=3e10, color=color.blue, make_trail=True, interval=10)
asteroid.v = vector(0,-0.75e4,0)
dt = 1e5


while True:
    rate(500)
    r = bigstar.pos - smallstar.pos
    F = G * bigstar.mass * smallstar.mass * r.hat / mag2(r)
    
    bigstar.p = bigstar.p - F*dt 
    smallstar.p = smallstar.p + F*dt

    bigstar.pos = bigstar.pos + (bigstar.p/bigstar.mass) * dt 
    smallstar.pos = smallstar.pos + (smallstar.p/smallstar.mass) * dt
    
    r1 = asteroid.pos - bigstar.pos
    r2 = asteroid.pos - smallstar.pos
    #print(r1)
    if (mag(r1)<3e9 or mag(r2)<3e9):
        print("collision")
        break
    elif (mag(r1)<8e9 or mag(r2)<8e9):
        dt = 4e3
        print("close call")
    elif (mag(r1)<8e10 or mag(r2)<8e10):
        dt = 2e4
    else: dt = 1e5
    asteroid.a = - G * bigstar.mass * r1.hat / mag2(r1) - G * smallstar.mass * r2.hat / mag2(r2)
    asteroid.v = asteroid.v + asteroid.a * dt
    asteroid.pos = asteroid.pos + asteroid.v *dt