from visual import *
from integrator import *
AU=149.6e9 #earth-sun orbital-radius
MU=384.4e6 #moon-earth orbital-radius

sun_mass=2e30
## keep sun motionless by inflating mass
sun_radius=6.96e8

earth_mass=6e24
earth_radius=6.37e6
earth_vel=2*math.pi*AU/(365.25*24.*60.*60.)

moon_mass=7.36e22
moon_radius=1.74e6
moon_vel=2*math.pi*MU/(27.32*24.*60.*60.)

boost=1.0
jupiter_mass=boost*1.9e27
jupiter_vel=2*math.pi*AU*5.2/(11.86*365.25*24.*60.*60)

scene.background=color.white
scene.autoscale=0
scene.range=6*AU

sun=sphere(pos=(0,0,0), velocity=vector(0,0,0),
             mass=sun_mass,radius=0.25*AU ,    color=color.yellow) 
earth=sphere(pos=(AU,0,0), velocity=vector(0,earth_vel,0),
             mass=earth_mass,radius=0.1*AU,   color=color.cyan)   
moon=sphere(pos=(AU+MU,0,0), velocity=vector(0,earth_vel+moon_vel,0),
             mass=moon_mass,radius=0.0005*AU,  color=color.red)
jupiter=sphere(pos=(5.2*AU,0,0),velocity=vector(0,jupiter_vel,0),
             mass=jupiter_mass, radius=0.15*AU, color=color.red)

# create list of gravitating objects
bodies=[sun,earth,jupiter]
for b in bodies:
    b.acc=vector(0,0,0)
    b.track=curve(color=b.color)

# set total momentum of system to zero
sum=vector(0,0,0)
for b in bodies:
    if (b!=sun):
        sum=sum+b.mass*b.velocity

sun.velocity=-sum/sun.mass

# dt corresponds to 3000 mins here
dt=30.*60.*100
while True:
    rate(100)
    update(bodies,dt)
    scene.center=sun.pos   #center the view on the sun