Почему я не могу вывести эту систему из двух N-тел без смещения (код части масс-центра не работает?)

Question

Я пытаюсь заставить эту диаграмму из 2 n тел работать в vpython, кажется, что она работает, но что-то не так с моим центром или массой, или что-то, я действительно не знаю. Система 2 n body смещается и не стоит на месте.

from visual import*
mt= 1.99e30 #Kg
G=6.67e-11 #N*(m/kg)^2

#Binary System stars

StarA=sphere(pos=  (3.73e11,0,0),radius=5.28e10,opacity=0.5,color=color.green,make_trail=True, interval=10)
StarB=sphere(pos=(-3.73e11,0,0),radius=4.86e10,opacity=0.5, color=color.blue,make_trail=True, interval=10)
#StarC=sphere(pos=(3.44e13,0,0),radius=2.92e11,opacity=0.5, color=color.yellow,make_trail=True, interval=10)

#mass of binary stars

StarA.mass= 1.45e30 #kg

StarB.mass= 1.37e30 #kg

#StarC.mass= 6.16e29 #kg

#initial velocities of binary system
StarA.velocity =vector(0,8181.2,0)
StarB.velocity =vector(0,-8181.2,0)
#StarC.velocity= vector(0,1289.4,0)

#Time step for each binary star

dt=1e5
StarA.pos= StarA.pos + StarA.velocity*dt
StarB.pos= StarB.pos + StarB.velocity*dt
#StarC.pos= StarC.pos + StarC.velocity*dt

#Lists

objects=[]
objects.append(StarA)
objects.append(StarB)
#objects.append(StarC)

#center of mass

Ycm=0
Xcm=0
Zcm=0
Vcmx=0
Vcmy=0
Vcmz=0
TotalMass=0
for each in objects:
   TotalMass=TotalMass+each.mass

   Xcm= Xcm + each.pos.x*each.mass
   Ycm= Ycm + each.pos.y*each.mass
   Zcm= Zcm + each.pos.z*each.mass

   Vcmx= Vcmx + each.velocity.x*each.mass
   Vcmy= Vcmy + each.velocity.y*each.mass
   Vcmz= Vcmz + each.velocity.z*each.mass

for each in objects:
    Xcm=Xcm/TotalMass
    Ycm=Ycm/TotalMass
    Zcm=Zcm/TotalMass
    Vcmx=Vcmx/TotalMass
    Vcmy=Vcmy/TotalMass
    Vcmz=Vcmz/TotalMass

each.pos=each.pos-vector(Xcm,Ycm,Zcm)
each.velocity=each.velocity-vector(Vcmx,Vcmy,Vcmz)

#Code for Triple system
firstStep=0   
while True:
rate(200)
for i in objects:
    i.acceleration = vector(0,0,0)
    for j in objects:
      if i != j:
        dist = j.pos - i.pos
        i.acceleration = i.acceleration + G * j.mass * dist / mag(dist)**3
for i in objects:
    if firstStep==0:
        i.velocity = i.velocity + i.acceleration*dt

        firstStep=1
    else:
        i.velocity = i.velocity + i.acceleration*dt


    i.pos = i.pos + i.velocity*dt

Answer 1

Два замечания, которые я хотел бы сделать по поводу вашей проблемы:

1. Вы используете метод интеграции Эйлера. Это часть вашего кода "i.velocity = i.velocity + i.acceleration * dt". Этот метод не очень точен, особенно с такими колебательными проблемами, как этот. Это одна из причин, по которой вы замечаете дрейф в вашей системе. Я бы рекомендовал использовать метод интеграции RK4. Это немного сложнее, но дает очень хорошие результаты, хотя интеграция с Verlet может быть более полезной для вас.

2. Ваша система может иметь чистый импульс в каком-то направлении. Они начинаются с одинаковой скорости, но имеют разные массы, поэтому суммарный импульс не равен нулю. Чтобы исправить это, преобразуйте систему в систему отсчета с нулевым импульсом. Для этого сложите импульс (вектор скорости и массы) всех звезд, а затем разделите их на общую массу системы. Это даст вам вектор скорости, который нужно вычесть из всех ваших начальных векторов скорости звезд, и ваш центр масс не будет двигаться.

Надеюсь это поможет!

Answer 2

Это не верно:

for each in objects:
    Xcm=Xcm/TotalMass
    ...

Это разделение должно произойти только один раз. И тогда средние должны быть удалены из объектов, как в

Xcm=Xcm/TotalMass
...
for each in objects:
    each.pos.x -= Xcm
    ...