title: Why I can't make this 2 N-body system orbit without shifting (Center of mass portion code not working?)
body: <p>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.   </p>
  
  <pre><code>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
  </code></pre>