#!/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