from __future__ import division
from visual import *

## Objects
Sun = sphere(pos=vector(0,0,0), radius=1e10, color=color.yellow)
Earth = sphere(pos=vector(1.5e11,0,0), radius=5e9, color=color.cyan)

# Constants
G = 6.67e-11
Sun.m = 2e30
Earth.m = 6e24
deltat = 20000

# Initial conditions
v = 2*pi*mag(Earth.pos)/(365*24*60*60)  #speed of the Earth in a circular orbit
Earth.p = Earth.m*vector(0, v, 0)
t = 0

#Visualization options
Earth.trail = curve(color=Earth.color)
scene.autoscale = 0

while t<10*(365*24*3600):
   rate(300)

   r = Earth.pos - Sun.pos
   rmag = sqrt(r.x**2 + r.y**2 + r.z**2)
   rhat=r/rmag

   Fmag = G*Earth.m*Sun.m/rmag**2
   Fnet = -Fmag*rhat

   Earth.p = Earth.p +Fnet*deltat

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

   t=t+deltat