# proj.py - External ballistics
# assumes the projectile has no pitch or yaw during the shot
# treats the projectile as a point mass--the CP is at the CG
# assumes forces other than drag and gravity are negligible
# assumes no spin of the projectile
# assumes no wind--any wind at all invalidates this analysis
from math import *

angle = raw_input("Angle (degrees) = ")
h = input("h (ft) = ")
Cd = input("Cd = ")
W = input("W (lb) = ")
Ain = input("A (in^2) = ")
v0 = input("v0 (ft/s) = ")

g = 32.174
rho = 0.0023283844
theta = pi * float(angle) / 180.0
m = W/g
A = Ain/144

def x1_dot():
	return -(0.5/(float(m)))*rho*Cd*A*x1*sqrt(pow(x1,2) + pow(y1,2))
	
def y1_dot():
	return -g - (0.5/(float(m)))*rho*Cd*A*y1*sqrt(pow(x1,2) + pow(y1,2))

def x0_dot():
	return x1

def y0_dot():
	return y1

x0 = 0
y0 = height
x1 = v0 * cos(theta)
y1 = v0 * sin(theta)

t = 0
# h = 0.01
h = 0.005

print "Calculating..."

ta = t
x0a = x0
y0a = y0
x1a = x1
y1a = y1

while y0a >= 0.0:
	# find x1a and y1a, the velocity components
	
	x1 = x1a
	y1 = y1a
	kx1a = x1_dot()
	ky1a = y1_dot()
	
	x1 = x1a + 0.5*h*kx1a
	y1 = y1a + 0.5*h*ky1a
	kx1b = x1_dot()
	ky1b = y1_dot()
	
	x1 = x1a + 0.5*h*kx1b
	y1 = y1a + 0.5*h*ky1b
	kx1c = x1_dot()
	ky1c = y1_dot()
	
	x1 = x1a + h*kx1c
	y1 = y1a + h*ky1c
	kx1d = x1_dot()
	ky1d = y1_dot()
	
	x1a = x1a + (float(h)/6)*(kx1a + 2*kx1b + 2*kx1c + kx1d)
	y1a = y1a + (float(h)/6)*(ky1a + 2*ky1b + 2*ky1c + ky1d)
	
	# find x0a and y0a, the position components
	
	x1 = x1a
	y1 = y1a
	kx0a = x0_dot()
	ky0a = y0_dot()
	
	x1 = x1a + 0.5*h*kx0a
	y1 = y1a + 0.5*h*ky0a
	kx0b = x0_dot()
	ky0b = y0_dot()
	
	x1 = x1a + 0.5*h*kx0b
	y1 = y1a + 0.5*h*ky0b
	kx0c = x0_dot()
	ky0c = y0_dot()
	
	x1 = x1a + h*kx0c
	y1 = y1a + h*ky0c
	kx0d = x0_dot()
	ky0d = y0_dot()
	
	x0a = x0a + (float(h)/6)*(kx0a + 2*kx0b + 2*kx0c + kx0d)
	y0a = y0a + (float(h)/6)*(ky0a + 2*ky0b + 2*ky0c + ky0d)
	
	# advance ta, the current time
	ta = ta + h

print "Done."

print "Hang time =", ta, "s"
print "Maximum range =", x0a, "ft"

vterm = sqrt(pow(x1a,2) + pow(y1a,2))

print "Terminal velocity =", vterm, "ft/s"

vave = float(x0a)/ta

print "Average speed =", vave, "ft/s"

mu = 3.82e-7
l = float(2)/12
Re = float(rho * l * vterm)/mu

print "Terminal Re =", Re

Re = float(rho * l * v0)/mu

print "Starting Re =", Re

ke = 0.5*m*pow(v0,2)

print "KEmax =", ke, "ft-lbf"

deltake = ke - 0.5*m*(pow(x1a,2) + pow(y1a,2))

print "deltaKE =", deltake,"ft-lbf"

losses = 100*float(deltake)/ke

print "% lost =", losses,"%"

ked = float(ke)/A

print "KEDmax =", ked, "ft-lbf/ft^2"

input("Exit...")