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

#/***************************************************************************
# *   Copyright (C) 2012, Paul Lutus                                        *
# *                                                                         *
# *   This program is free software; you can redistribute it and/or modify  *
# *   it under the terms of the GNU General Public License as published by  *
# *   the Free Software Foundation; either version 2 of the License, or     *
# *   (at your option) any later version.                                   *
# *                                                                         *
# *   This program is distributed in the hope that it will be useful,       *
# *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
# *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
# *   GNU General Public License for more details.                          *
# *                                                                         *
# *   You should have received a copy of the GNU General Public License     *
# *   along with this program; if not, write to the                         *
# *   Free Software Foundation, Inc.,                                       *
# *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
# ***************************************************************************/
 
import sys, math

# universal gravitational constant
G = 6.67428e-11

# earth mass
e_mass = 5.9736e+024

orbit_radius = 8.0e7

# circular orbit velocity
cov = math.sqrt(G * e_mass / orbit_radius)

# position
pos = complex(orbit_radius,0)

# velocity: set elliptical orbit: cov/3
vel = complex(0,cov)/3

dt = .1 # delta-t seconds

n = 0

# energy max and min
emax = -1e9
emin = 1e9

# simulate orbit

while True:
  # update velocity
  vel += pos * -G * e_mass * dt * abs(pos)**-3
  # update position
  pos += vel * dt
  # get magnitudes
  r = abs(pos)
  v = abs(vel)
  # calculate total energy: KE + PE
  e = (v**2)/2 -(G * e_mass / r)
  # collect min and max energy bounds
  emax = max(emax,e)
  emin = min(emin,e)
  # compute error
  error = abs((emax-emin) * 100 / emin)
  # display results
  # show only a few results
  if(n % 10000 == 0):
    print 'p = %8.4e v = %8.4e e = %8.4e error = %6.2e%%' % (r,v,e,error)
  n += 1