from visual import *


G               = 6.67428e-11       # N m^2 / kg^2
sun_mass        = 1.98892e30        # kilograms

# Properties of the Earth's orbit
eccentricity    = 0.016710219
a               = 1.495978875e11    # meters
year            = 3.1556926e7       # seconds

pos_perihelion = (1-eccentricity)*a
vel_perihelion = sqrt(((1+eccentricity)*G*sun_mass)/((1-eccentricity)*a)) # m/s

# Define a window variable so we can shut off autoscale later.
window = display(title="Sun and Earth", autoscale=1)


# Create the Sun.
sun = sphere()
sun.pos = vector(0,0,0)
sun.radius = 50*6.955e8 # meters
sun.color = color.yellow
sun.mass = sun_mass
sun.p = vector(0, 0, 0) * sun.mass

# Create the Earth.
earth = sphere()
earth.pos = vector(pos_perihelion,0,0)
earth.radius = 2000*6.4e6
earth.color = color.blue
earth.mass = 5.9736e24
earth.p = vector(0,vel_perihelion, 0) * earth.mass

# Shut off autoscale so the screen doesn't flickr.
window.autoscale = 0

for a in [sun, earth]:
  a.orbit = curve(color=a.color, radius = earth.radius/10)

time    = 0
dt      = 60.0 # seconds

dist    = earth.pos - sun.pos
maxdist = abs(dist)
mindist = abs(dist)

while time < year:
    
  # rate(100)

  dist = earth.pos - sun.pos
  if abs(dist) > maxdist:
      maxdist = abs(dist)
      
  force = G * sun.mass * earth.mass * dist / mag(dist)**3
  sun.p = sun.p + force*dt
  earth.p = earth.p - force*dt

  time = time + dt

  print time

  for a in [sun, earth]:
    a.pos = a.pos + a.p/a.mass * dt
    a.orbit.append(pos=a.pos)

    
print "Aphelion:", maxdist