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