from visual import *


G               = 6.6742867e-11        # N m^2 / kg^2 (http://en.wikipedia.org/wiki/Gravitational_constant)
#sun_mass        = 1.98892e30        # kilograms (http://en.wikipedia.org/wiki/Solar_mass)
sun_mass        = 1.989100e30        # kilograms (http://nssdc.gsfc.nasa.gov/planetary/factsheet/sunfact.html)


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

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

r_initial       = 0.983 * a         # meters
elat_initial    = 0.00              # degrees
elon_initial    = 100.68            # degrees
vel_initial     = vel_perihelion    # THIS NEEDS TO BE CORRECTED!!!!

x_initial       = r_initial * cos(elon_initial*pi/180) * cos(elat_initial*pi/180)
y_initial       = r_initial * sin(elon_initial*pi/180) * cos(elat_initial*pi/180)
z_initial       = r_initial * sin(elat_initial*pi/180)


# 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(x_initial,y_initial,z_initial)
earth.radius = 2000*6.4e6
earth.color = color.blue
earth.mass = 5.9736e24
earth.p = vector(0,vel_initial, 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      = 86400.0 # seconds

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

while time < year:

  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