from visual import *

#Constants
G = 6.7e-11
day = 24*60*60
year = 365*day
jyear = 11.86*year
rscale = 3
mscale = 200

#Objects
Sun = sphere(pos=vector(0,0,0), radius=4e10, color=color.yellow)
Jupiter = sphere(pos=vector(8e11,0,0), radius=3e10, color=color.orange)
Earth = sphere(pos=vector(rscale*1.5e11,0,0), radius = 2e10, color=color.cyan)
#scene.autoscale = False

Sun.m = 2e30
Jupiter.m = 2e27*mscale
Earth.m = 6e24

#Initial values
Earth.p = Earth.m*vector(0,2*pi*Earth.pos.x/year/rscale**1.5,0)
Earth.trail=curve(color=Earth.color)
Earth.trail.append(pos=Earth.pos)

Jupiter.p = Jupiter.m*vector(0,2*pi*Jupiter.pos.x/jyear,0)

Sun.p= -(Earth.p+Jupiter.p)

deltat = day
t=0

#Loop for repetitive calculations
while t < 100*year:
    rate(2000)
    rES = Earth.pos - Sun.pos
    rESmag = mag(rES)
    rEShat = rES/rESmag
    FnetES = (-G*Earth.m*Sun.m/rESmag**2)*rEShat

    rEJ = Earth.pos - Jupiter.pos
    rEJmag = mag(rEJ)
    rEJhat = rEJ/rEJmag
    FnetEJ = (-G*Earth.m*Jupiter.m/rEJmag**2)*rEJhat

    rJS = Jupiter.pos - Sun.pos
    rJSmag = mag(rJS)
    rJShat = rJS/rJSmag
    FnetJS = (-G*Jupiter.m*Sun.m/rJSmag**2)*rJShat

    FnetEarth = FnetES + FnetEJ
    FnetJupiter = FnetJS - FnetEJ
    FnetSun = -(FnetES + FnetJS)

    Earth.p = Earth.p + FnetEarth*deltat
    Earth.pos = Earth.pos+Earth.p/Earth.m*deltat
    Earth.trail.append(pos=Earth.pos)

    Jupiter.p = Jupiter.p + FnetJupiter*deltat
    Jupiter.pos = Jupiter.pos + Jupiter.p/Jupiter.m*deltat

    Sun.p = Sun.p + FnetSun*deltat
    Sun.pos = Sun.pos + Sun.p/Sun.m*deltat

    t = t + deltat