我这里有一段代码,可以发现两个粒子何时彼此近距离接触。 5个粒子呢?

I have here a block of code that finds when two particles have a close encounter with each other. What about 5 particles?

这个代码的原因是为了找到一个行星系统有2个行星(粒子)彼此近距离相遇所需的时间。但是,我想解决具有 5 个粒子的行星系统的时间问题。这将打印在最后一行 print(sim.t).

http://rebound.readthedocs.io/en/latest/examples.html 这里有一些关于 REBOUND 用法的文档,您可以在其中找到对我所做的一些声明的最佳解释。

sim = setupSimulation()
sim.exit_min_distance = 0.01
Noutputs = 10000
year = 2.*np.pi
times = np.linspace(0.,10E+9.*year, Noutputs)
distances = np.zeros(Noutputs)
ps = sim.particles
try:
    for i,time in enumerate(times):
        sim.integrate(time)
        dp = ps[1] - ps[2]   
        distances[i] = np.sqrt(dp.x*dp.x+dp.y*dp.y+dp.z*dp.z)
except rebound.Encounter as error:
    print(error)

print(sim.t)

第 2 行定义了两个粒子之间的近距离接触(达到此值时模拟将停止)。

第11行取粒子坐标差,如果<=.01,则打印近距离接触

我在考虑 if-then 语句:

if ps[1] - ps[2] <= .01:
    dp = ps[1] - ps[2]
else ps[2] - ps[3] <= .01:
    dp = ps[2] - ps[3]

等等...

我想在我 运行 之前确保它能正常工作,因为模拟时间是 10e+9 年,并且只能 运行 在本地超级计算机上获得结果合理的时间。

我不明白为什么 if 语句方法无法检测任何两个粒子是否彼此足够接近。编辑:编写此示例就像您的 2 个粒子代码示例正在运行一样。

根据您是否需要检查非连续粒子,您需要生成每个组合并进行检查(如示例所示)。

看来您还需要不同的距离枚举策略。如果索引不必与时间计数器匹配,您可以增加一个单独的计数器(如示例所示)。

您还可以使用 itertools 生成 if 语句以避免拼写错误并提高清晰度:

import itertools

# other setup code

counter = 0
try:
    for time in times:
        sim.integrate(time)
        for p0, p1 in itertools.combinations(range(1, 6)):  # generates combinations of the particles numbered 1 thru 5
            dp = ps[p0] - ps[p1]
            distances[counter] = np.sqrt(dp.x*dp.x+dp.y*dp.y+dp.z*dp.z)
            counter += 1
except rebound.Encounter as error:
    print(error)

我再补充一点——为了确保您的代码有效,您可能应该为它编写某种测试套件。如果它解决了与从不解决会发生什么。

您可以创建一个对象来模拟 sim 及其方法,在您 知道 将满足或不满足退出条件的位置返回主体,而无需所有开销需要一台超级计算机。您可以将它们硬编码到一个小列表中。