#!/usr/bin/env python
# -*- coding: utf-8 -*-

import math as m

import matplotlib.pyplot as plt

plt.axes().set_aspect('equal')
plt.xlabel("X values")
plt.ylabel("Y values")
plt.title("Orbit Example")
plt.grid()

# "Big G", universal gravitational constant
# https://physics.nist.gov/cgi-bin/cuu/Value?bg

G = 6.67430e-11

# Sun: https://nssdc.gsfc.nasa.gov/planetary/factsheet/sunfact.html

sunmass = 1.9885e30  # kilograms

# Earth: https://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html

earth_aphelion = 149.6e9  # greatest distance from the sun, meters
earth_velocity = 29.78e3  # initial velocity m/s

# model time slice seconds, default one hour

Δt = 3600.0

# x and y initial positions

yp = oy = 0
xp = earth_aphelion

# x and y initial velocities

yv = -earth_velocity
xv = 0

# arrays to contain orbital data

xvals = []
yvals = []

# precompute an acceleration term

accel_term = sunmass * -G * Δt

# model the orbit

while oy <= 0 or yp > 0:
  oy = yp
  
  # acquire gravitational acceleration
  acc = m.pow(xp**2+yp**2, -1.5) * accel_term
  
  # update velocity
  xv += xp * acc
  yv += yp * acc
  
  # update position
  xp += xv * Δt
  yp += yv * Δt

  # save values for later plotting
  xvals += [xp]
  yvals += [yp]

# draw the sun
plt.scatter(0, 0, marker='o', s=240, color='#ffff20', edgecolor='black')

# plot the orbit
plt.plot(xvals, yvals, color='g')

plt.savefig("gravity_example1.png")