# bags.py - btrettel's air gun simulation
from math import *

# inputs
Co = 6 # flow coefficient at no pressure differential (typically called Cv)
d = 1.590 # barrel diameter, in
L = 18 # barrel length, in
P = 50 # air chamber pressure, psig
Pf = 2 # minimum pressure to move projectile, psia
Tc = 530 # initial air chamber temperature, R
Vc = 40 # air chamber volume, in^3
Vd = 10 # dead volume, in^3
Vp = 0.1 # pilot volume, in^3
W = 0.0749571691 # weight of projectile, lb
to = 5e-3 # valve opening time, s
n = 0 # fudge factor to account for changes in flow coefficient with pressure
      # a constant Cv has an n of 0

# constants
Cp = 2241.1 # specific heat at constant pressure, in/R
Cv = 1606.1 # specific heat at constant volume, in/R
g = 32.2 # gravitational acceleration, ft/s^2
k = 1.4 # ratio of the specific heats
lbtocg = 45360 # multiply by something in pounds to convert to centigrams
Patm = 14.7 # atmospheric pressure, psia
R = 640.12 # gas constant for air, in/R
Tr = 528 # SCFM reference temperature, R
T = 530.0 # atmospheric temperature, R
SG = 1 # specific gravity of the gas at atmospheric pressure and 60 degrees F
pr = 0.528 # critical pressure ratio

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

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

A = (pi/4) * pow(d,2)

# functions
def C(Pi):
	return Co * pow(Pi,n)

def Vb(xi):
	# finds the current barrel volume
	return (pi/4) * pow(d,2) * xi + Vd

def Q(Pi, Po, Tci):
	# finds the volumetric flow rate at STP
	if pr * Pi <= Po:
		flow = 462.24 * C(Pi - Po) * sqrt((pow(Pi,2) - pow(Po,2))/(Tci * SG))
	else:
		flow = (316.8 * C(Pi - Po) * Pi)/sqrt(SG * Tci)
	
	if t < to:
		return flow * (t / to)
	else:
		return flow

def x_2(Pbi,Tbi):
	# finds the current projectile acceleration
	if x_dot > 0:
		Pfi = Pf
	else:
		Pfi = -Pf
	force = ((12 * pi * pow(d,2))/4) * (Pbi + ((0.5 * 12)/g) * (Pbi/(R * Tbi)) * pow(x_dot/12,2) - Patm - ((0.5 * 12)/g) * (Patm/(R * T)) * pow(x_dot/12,2))
	friction = ((12 * pi * pow(d,2))/4) * Pf
	if abs(force) > abs(friction):
		return (force - friction)/m
	else:
		return 0

def Tt(mi,mo,Ti,To,Poi,worki):
	# finds the new temperature with mixing, assuming no heat transfer to the surroundings
	return (mi * Ti * h + mo * To)/(mi * h + mo) - (worki)/(Cv * (mi * h + mo))

def m_dot(Qi):
	# find mass flow rate from Q
	return (Patm * Qi)/(R * Tr)

def Ps(mi, Ti, Vi):
	# find current pressure
	return (mi * R * Ti)/Vi

def E(Pmax,Pmin,Vmin):
	# figures out energy stored in a gas
	return (1.0/12) * (Pmin * pow(Pmax/Pmin,1/k) * Vmin - Pmax * Vmin)/(1 - k)

def loop():
	global Pb, Pc, Tb, Tc, mb, mc, h, x_dot, x, work, t
	Pbo = Pb
	Pco = Pc
	Tbo = Tb
	Tco = Tc
	
	if Pb <= Pc: # if barrel pressure < chamber pressure --> flow into the barrel
		m_dota = m_dot(Q(Pc, Pb, Tc))
		Tb = Tt(m_dota,mb,Tc,Tb,Pb,work)
		mb = mb + m_dota * h
		mc = mc - m_dota *h
	else: # if barrel pressure > chamber pressure --> flow into the chamber
		m_dota = m_dot(Q(Pb, Pc, Tb))
		Tc = Tt(m_dota,mc,Tb,Tc,Pc,0)
		mb = mb - m_dota * h
		mc = mc + m_dota *h
	
	Pb = Ps(mb, Tb, Vb(x))
	Pc = Ps(mc, Tc, Vc)
	
	x_dot_0 = x_dot
	
	accel = x_2(Pb,Tb)
	x_dot = x_dot + h * accel
	x = x + h * x_dot
	
	work = (0.5/12) * m * (pow(x_dot,2) - pow(x_dot_0,2))
	
	t = t + h

# initial conditions before looping
t = 0
x = 0
work = 0
x_dot = 0
Pb = Patm
Pc = P + Patm
Tb = T
mb = (Pb * Vb(x))/(R * Tb)
mc = (Pc * Vc)/(R * Tc)

while x < L and t < (0.1 - h):
	loop()

print "Vf =", x_dot/12
print "tf =", t
mair = lbtocg * (P * (Vc + Vp))/(R * T)
print "Moved air mass, with pilot volume (cg) =", mair
mair3 = lbtocg * (P * Vc)/(R * T)
print "Moved air mass, zero pilot volume (cg) =", mair3
mair2 = lbtocg * ((P - Pc + Patm) * Vc)/(R * T)
print "Moved air mass, new valve (cg) =", mair2
print "Old valve vs. no pilot  =", (1 - mair3/mair)
print "New valve vs. old  =", (1 - mair2/mair)
print "New valve vs. no pilot old  =", (1 - mair2/mair3)
ke = 0.5 * pow(x_dot/12,2) * m
print "KE =", ke
pe = E(P + Patm,Patm,Vc + Vp)
print "PE =", pe
pes = E(Pc,Patm,Vc)
print "PE saved =", pes
print "eta =", ke/pe
pe2 = E(P + Patm,Patm,Vc)
print "eta (no pilot) =", ke/pe2
print "eta2 =", ke/(pe2 - pes)
print "	Barrel		Chamber"
print "T:" + "\t" + str(Tb-459.67) + "\t" + str(Tc-459.67)
print "P:" + "\t" + str(Pb-Patm) + "\t" + str(Pc-Patm)