title: Python gravity simulator behaving strangely
body: <p>I'm making a gravity simulation in Python (in 3D with VPython, to be exact) and I'm sure there's nothing wrong with the code, but it behaves strangely when two objects get close to each other.</p>
  
  <p>My inspiration is <a href="http://testtubegames.com/gravity.html" rel="nofollow">http://testtubegames.com/gravity.html</a>. Note how you can place two planets with no velocity, they move towards each other, overtake, decelerate and turn back. In my program, they overtake, and decelerate, but only proportionately to the distance, so technically it should never turn back anyway.</p>
  
  <p>I realize that the law f=G*(m1*m2)/r**2 wouldn't work if r (the distance) is or gets too close to 0, so I have included a max-out, so if it is less than 1 it is set to 1 (units not in pixels, by the way), but it still does not work.</p>
  
  <p>Simple logic also suggests that the objects should not react in this way, so the next thing that follows is that I must be missing something.</p>
  
  <p>Here is an extract of the code:</p>
  
  <pre><code>from visual import *
  a = sphere(x=-10,mass=10, vel=vector())
  b = sphere(x=10, mass=10, vel=vector())
  
  while 1:
      rate(20)
  
      #distance between the two objects, a and b, where a.r.mag would be the magnitude of the vector
      a.r = b.pos - a.pos
      b.r = a.pos - b.pos
  
      a.force = a.r
      if a.r.mag &gt; 1:
          a.force.mag = (a.mass * b.mass) / a.r.mag**2
      else:
          a.force.mag = (a.mass * b.mass) / 1
      a.vel = a.vel + a.force / a.mass
  
  
      b.force = b.r
      if b.r.mag &gt; 1:
          b.force.mag = (a.mass * b.mass) / b.r.mag**2
      else:
          b.force.mag = (a.mass * b.mass) / 1
      b.vel = b.vel + b.force / b.mass
  
      a.pos = a.pos + a.vel
      b.pos = b.pos + b.vel
  </code></pre>
  
  <p>EDIT: Code re-written in response to shockburner:</p>
  
  <pre><code>from visual import *
  import sys
  
  limit2 = sys.float_info.min
  limit = limit2**0.5
  timestep = 0.0005
  
  a = sphere(x=-5,mass=10, vel=vector())
  b = sphere(x=5, mass=10, vel=vector())
  
  def force(ob1, ob2):
      ob1.r = ob2.pos - ob1.pos
      ob1.force = ob1.r + vector()
      if ob1.r.mag &gt; limit:
          ob1.force.mag = (ob1.mass * ob2.mass) / ob1.r.mag2
      else:
          ob1.force.mag = (ob1.mass * ob2.mass) / limit2
      return ob1.force
  
  while 1:
      rt = int(1/timestep)
      rate(rt)
  
      a.acc = force(a, b) / a.mass
      b.acc = force(b, a) / b.mass
  
      a.pos = a.pos + timestep * (a.vel + timestep * a.acc / 2)
      b.pos = b.pos + timestep * (b.vel + timestep * b.acc / 2)
  
      a.acc1 = force(a,b) / a.mass
      b.acc1 = force(b,a) / b.mass
  
      a.vel = a.vel + timestep * (a.acc + a.acc1) / 2
      b.vel = b.vel + timestep * (b.acc + b.acc1) / 2
  </code></pre>
  
  <p>Any help or pointer in the right direction would be greatly appreciated, and if the answer turns out to be idiotically simple (which with me is usually the case) remember I am quite an idiot anyway.</p>