# bags_0.1.py - Internal ballistics of a pneumatic gun
# 12/29/08
from math import *

# inputs

# Vc = input("Volume of air chamber (in^3) = ")
# L = input("Barrel length (in) = ")
# d = input("Barrel diameter (in) = ")
# P = input("Pressure (psig) = ")
# W = input("Projectile weight (lb) = ")
# Pf = input("Equivalent pressure of friction (psi) = ")
# Vd = input("Dead volume (in^3) = ")
# Cv = input("Cv = ")
# Vp = input("Pilot volume (in^3) =")

Vc = 40
L = 18
d = 1.590
P = 50
W = 0.0749571691
Pf = 2
Vd = 10
Cv = 6
Vp = 0.1

# program specs
h = 1e-4 # interval in seconds

# constants
T = 530 # atmospheric temperature in degrees Rankine
Patm = 14.7 # atmospheric pressure in psia
g = 32.2 # acceleration due to gravity in ft/s^2

# convert P in gauge pressure to absolute pressure
Pc = Patm + P;

# convert weight W in lb to mass m in slug
m = W/g;

def Vb():
	# Finds current barrel volume from X and d
	return (pi/4) * pow(d,2) * x + Vd

# before looping...
t = 0
x = 0
x_dot = 0
Pb = Patm
PVb = Pb * Vb()
PVc = Pc * Vc

while x < L and t < (0.1 - h):
	# keep accelerating the projectile until it leaves the barrel
	if Pc < Pb:
		#Pc = (Pb*Vb() + Pc*Vc) / (2*(Vb()+Vc))
		Pb = Pc
	PVb_dot = 462.24 * Cv * sqrt((pow(Pc,2) - pow(Pb,2))/T) * Patm
	PVb = PVb + PVb_dot * h
	PVc = PVc - PVb_dot * h
	Pb = PVb/Vb()
	Pc = PVc/Vc
	x1_dot = ((12 * pi * pow(d,2))/(m * 4)) * (Pb - Patm - Pf)
	x_dot = x_dot + h * x1_dot
	x = x + h * x_dot
	t = t + h

R = 640.12 # gas constant for air, in/R
lbtocg = 45360 # multiply by something in pounds to convert to centrigrams
print "Vf =", x_dot/12
print "Tf =", t
print "Pf =", Pb
ke = 0.5 * pow(x_dot/12,2) * m
pe = (1.0/12) * (P + Patm) * (Vc + Vp) * log((P + Patm)/Patm)
print "KE =", ke
print "PE =", pe
print "eta =", ke/pe
input("Done.")