from visual import *

gravG=6.7e-11 #gravitational const in MKS units, which are used throughout

def Eulerstep0(y,h,f):
#This version does not use lists.
    return y +f*h

#Initialize giant's properties: position, size, color, mass, momentum
giant = sphere()
giant.pos = vector(-1e11,0,0)
giant.radius = 2e10
giant.color = color.red
giant.mass = 2e30
giant.p = vector(0, 0, -1e4) * giant.mass #Here p is momentum.  Note 
#that initially moving perpendicular to vector between stars.

#Similarly, initialize dwarf
dwarf = sphere()
dwarf.pos = vector(1.5e11,0,0)
dwarf.radius = 1e10
dwarf.color = color.yellow
dwarf.mass = 1e30
dwarf.p = -giant.p

for astar in [giant, dwarf]:  #Note use of list for iteration
  astar.orbit = curve(color=astar.color, radius = 2e9)
#set up plotting of orbit with VPython curve

dt = 86400

while 1:
  rate(100)

  dist = dwarf.pos - giant.pos #points giant to dwarf
  force = -gravG * giant.mass * dwarf.mass * dist / mag(dist)**3 #attractive
  dwarf.p = Eulerstep0(dwarf.p,dt,+force)   #Euler stepping for momentum
  giant.p = Eulerstep0(giant.p,dt,-force)   #equal and opposite

  for astar in [giant, dwarf]:
    astar.pos = astar.pos + astar.p/astar.mass * dt  #Euler step for position
    astar.orbit.append(pos=astar.pos)  #add to plotted curve