# Remove if not using Blender	
#import Blender
#from Blender import *

import math
from math import pi as pi
from math import *

# Remove if not using Blender	
#scn = Blender.Scene.GetCurrent()


# Let 1 Blender unit = 10 km
# Work in SI units but then convert for diaply

Mearth = 5.972E24
G = 6.67E-11
Rearth = 6371E3
BU = 1.0 / 10E3
km = 1000.0
hour = 3600.0
day = 86400.0
week = 604800.0
year = 31449600.0
gyr  = 3.15569E16
kpc  = 3.08567758E19

# USER INPUTS
# Initial position in km
x = 0.0
y = 30000.0*km	#Rearth*2
# Initial Velocity in km/s
vi = -2.181
# Direction of velocity (angle from vertical in degrees)
vtheta = 90.0
# Timestep in seconds. Don't go below 0.5 or numerical errors dominate.
tstep = 1.0
# Maximum time in seconds
tmax = 7*hour
# Maximum number of timesteps
nmax = 40000
# If using Blender, probably not a good idea to animate every single timestep
# - this will generate a huge, slow file. Better to only output every specified
# number of steps.
# The Earth rotation is currently animated (somewhat imprecisely) assuming a 1 
#second timestep outputevery 51 frames.
fstep = 51

# END OF USER INPUTS


# Initial radial position. Will stop when within the radius of Earth
rpos = sqrt(x**2 + y**2)

# Resolve the initial velocity into x and y components, convert to SI
vtheta = radians(90.0-vtheta)
vx = vi*cos(vtheta) * 1000.0
vy = vi*sin(vtheta) * 1000.0

# Animate this object. Remove if not using Blender	
Sat = Blender.Object.Get('Satellite5')


# Calculate x and y accleration given current position
def acc(xpos, ypos):
	r = sqrt(xpos**2 + ypos**2)
	a = (G*Mearth) / r**2
	
	if x <> 0.0:
		theta = atan(abs(ypos)/abs(xpos))
	if x == 0:
		theta = pi / 2.0
	
	# Accelerations will be negligible if x==0 or y==0, so no point adding
	# checks for these cases.	
	if (xpos >= 0.0):
		ax = -a*cos(theta)
	if (xpos < 0.0):
		ax = a*cos(theta)
	if (ypos >=0.0):
		ay = -a*sin(theta)
	if (ypos < 0.0):
		ay = a*sin(theta)
		
	return ax, ay


# Remove if not using Blender	
me = Blender.Mesh.New('Track')


# Count the number of steps
n = 0

# Simulation time in seconds
t = 0

frame = 1

while (((rpos) > Rearth) and (n < nmax) and (t < tmax)):
	# Calculate the distance travelled, assuming constant acceleration for
	# the specified timestep
	# First calculate the initial acceleration
	
	ax, ay = acc(x,y)
		
	x = x + vx + (0.5*ax*(tstep**2))
	y = y + vy + (0.5*ay*(tstep**2))
	
	vx = vx + (ax*tstep)
	vy = vy + (ay*tstep)
	
	# Generate the track. Remove if not using Blender	
	me.verts.extend(x*BU,y*BU,0.0)
	
	
	n = n + 1
	
	rpos = sqrt(x**2 + y**2)
	
	t = t + tstep
	
	# Animates the Blender object. Remove or comment out if not using Blender.
	if (int(t)==1) or (t % 51 == 0):
		frame = frame + 1
		Blender.Set('curframe',frame)
	
		Sat.LocX = x*BU
		Sat.LocY = y*BU
	
		Sat.insertIpoKey(Blender.Object.LOC)
	
	print t/3600.0

# Remove if not using Blender	
ob = scn.objects.new(me,'Track')
	

print 'Time elapsed : ', t/hour, 'hours.'
