Python重力模拟器表现怪异

我正在用Python(确切地说是用VPython的3D)进行重力模拟,我确信代码没有错,但当两个物体靠近时,它的行为会很奇怪。

我的灵感来自 http://testtubegames.com/gravity.html 。注意如何放置两颗没有速度的行星,它们朝着对方移动,超越,减速并返回。在我的程序中,它们会超车和减速,但只与距离成比例,所以从技术上讲,无论如何都不应该回头。

我意识到,如果r(距离)等于或过于接近0,则f=G*(m1*m2)/r**2定律将不起作用,所以我包含了一个max-out,所以如果它小于1,则将其设置为1(顺便说一句,单位不是像素),但它仍然不起作用。

简单的逻辑也表明对象不应该以这种方式做出反应,所以接下来的事情是我一定错过了什么。

以下是代码摘录:

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 > 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 > 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

编辑:针对冲击燃烧器重新编写的代码:

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 > 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

任何正确方向的帮助或指示都将不胜感激,如果答案简单得愚蠢(对我来说通常是这样),请记住,无论如何,我都是个白痴。