Why I can't make this 2 N-body system orbit without shifting (Center of mass portion code not working?)

Question

I am trying to make this 2 n body diagram to work in vpython, it seems that is working but something is wrong with my center or mass, or something, i don't really know. The 2 n body system is shifting and is not staying still.

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

Two points I'd like to make about your problem:

1. You're using Euler's method of integration. That is the "i.velocity = i.velocity + i.acceleration*dt" part of your code. This method is not very accurate, especially with oscillatory problems like this one. That's part of the reason you're noticing the drift in your system. I'd recommend using the RK4 method of integration. It's a bit more complicated but gives very good results, although Verlet integration may be more useful for you here.

2. Your system may have a net momentum in some direction. You have them starting with the same speed, but have different masses, so the net momentum is non-zero. To correct for this transform the system into a zero momentum reference frame. To do this add up the momentum (mass times velocity vector) of all the stars, then divide by the total mass of the system. This will give you a velocity vector to subtract off from all your initial velocity vectors of the stars, and so your center of mass will not move.

Hope this helps!

Answer 2

This is wrong:

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

This division should only happen once. And then the averages should be removed from the objects, as in

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