# BAGS - http://trettel.org/bags/
# btrettel's air gun simulation
# bags.py - standalone script
# Last updated 07/16/09
# version 1.0.1b

# BAGS assumes a perfect gas (constant specific heats plus ideal gas)
# BAGS neglects heat transfer and shock formation
# BAGS is inaccurate above Mach 0.5 or so

# A good place to learn about compressible flow from is NASA
# See http://www.grc.nasa.gov/WWW/K-12/airplane/shortc.html
# and http://www.grc.nasa.gov/WWW/BGH/shorth.html
# Most of my links in the code are to Wikipedia, but NASA is very good too

# BAGS requires the math and csv Python modules

from math import *
import csv

# inputs that describe the gun

gas = 'air'
Cva = 2.5 # flow coefficient (typically called Cv)
d = 13.4e-3 # barrel diameter, m
L = 0.3048 # barrel length, m
m = 1.3e-3 # mass of projectile, kg
P = 5e5 # gauge pressure, Pa
Pf = 14e3 # minimum pressure to move projectile, Pa
to = 5e-3 # valve opening time, s
Vc = 3.3e-5 # gas chamber volume, m^3
Vd = 8.2e-6 # dead volume, m^3
Vp = 1.6e-6 # pilot volume, m^3

# configuration

Pr = 101325.0 # reference/standard atmospheric pressure in Pa
Rbar = 8314.472 # gas constant in J/kmol*K
g = 9.80665 # gravitational acceleration in m/s^2
Tr = 293.15 # reference temperature in K
h = 1e-6 # time interval in seconds

# atmospheric conditions (pressure and temperature) available from
# http://weather.gov in the USA and a barometer there and elsewhere; see
# http://en.wikipedia.org/wiki/Inch_of_mercury to convert barometer measurements
# use Tatm = Tr for the reference temperature (293.15 K / 20.15 C / 68 F)
# use Patm = Pr for the standard atmospheric pressure
# use Tc = Tatm and Tb = Tatm unless gas in gun is hot or cold

Tatm = 293.0 # atmospheric temperature in K
Tc = Tatm # initial gas chamber temperature, K
Tb = Tatm # initial barrel gas temperature, K
Patm = Pr # atmospheric temperature in Pa

def gasprop(prop,gas):
	# reads gas properties from gas.txt
	# syntax is gasprop('desired property','gas name')
	# ex. gasprop('k','N2')
	
	reader = csv.reader(open("gas.txt", "rb"), delimiter = '	')
	
	data = ['Cp', 'Cv', 'MW']
	
	item = data.index(prop) + 1
	
	for row in reader:
		if row[0][0] != '#':
			if row[0] == gas:
				return float(row[item])

def prcrit(k):
	# returns the critical pressure ratio, above which flow velocity is
	# limited to the sonic velocity
	# k is the ratio of specific heats of the gas
	# See http://en.wikipedia.org/wiki/Choked_flow for discussion
	
	return pow(2 / (k + 1), k / (k - 1))

def a(k, MW, T, Rbar):
	# returns the sonic velocity in m/s
	# k is the ratio of specific heats of the gas
	# MW is the molecular weight of the gas in kg/kmol
	# T is the temperature of the gas in K
	# Rbar is the gas constant in J/kmol*K
	# Is the simple speed of sound equation for a perfect gas
	# See http://en.wikipedia.org/wiki/Speed_of_sound#Speed_in_ideal_gases_and_in_air
	
	return sqrt(k * Rbar * T / MW)

def Pt(P, k, v, MW, T, Rbar):
	# returns the total pressure (static and dynamic added) in Pa
	# P is the static pressure in Pa
	# k is the ratio of specific heats of the gas
	# v is the velocity of the gas in m/s
	# MW is the molecular weight of the gas in kg/kmol
	# T is the temperature of the gas in K
	# Rbar is the gas constant in J/kmol*K
	# Is the simple total pressure equation for a compressible fluid
	# See http://en.wikipedia.org/wiki/Stagnation_pressure#Compressible_flow
	
	return P * pow(1 + (pow(v, 2) * MW * (k - 1)) / (2 * Rbar * k * T), k / (k - 1))

def m_dot(Pin, Pout, Tin, k, t, to, Cva, SG, MW, pr):
	# returns the mass flow through the valve in kg/s
	# Pin is the pressure of the gas stream flowing in in Pa
	# Pout is the pressure of the gas stream flowing out in Pa
	# Tin is the temperature of the gas stream flowing in in K
	# k is the ratio of specific heats of the gas
	# t is the current time in s
	# to is the opening time of the valve in s
	# Cva is the valve flow coefficient
	# SG is the specific gravity of the gas
	# MW is the molecular weight of the gas in kg/kmol
	# pr is the critical pressure ratio
	
	# Determine if the flow is sonic (which occurs if the pressure ratio
	# is above the critical pressure ratio) and from there calculate the
	# flow
	# I took these equations from a book and do not know what assumptions
	# they make
	# Similar equations can be found at http://www.engineeringtoolbox.com/flow-coefficients-d_277.html
	# In the future I plan to make this part more robust by modeling it as
	# a hole in a plate and putting a force balance on the piston
	
	if pr * Pin <= Pout:
		flow = 8.189e-7 * Cva * sqrt((pow(Pin, 2) - pow(Pout, 2)) / (Tin * SG))
	else:
		flow = (5.612e-7 * Cva * Pin) / sqrt(SG * Tin)
	
	# if the valve is still opening, scale the flow linearly against the
	# fraction the valve is open
	
	if t < to:
		flow = flow * (t / to)
	
	return (Patm * flow * MW) / (Rbar * Tr)

def frange(start, end=None, inc=None):
    # a range function that accepts float increments

    if end == None:
        end = start + 0.0
        start = 0.0

    if inc == None:
        inc = 1.0

    L = []
    while 1:
        next = start + len(L) * inc
        if inc > 0 and next >= end:
            break
        elif inc < 0 and next <= end:
            break
        L.append(next)
        
    return L

def Vb(x, A, Vd):
	# returns the barrel volume at projectile position x (m) in m^3
	# A is the cross-sectional area of the barrel in m^2
	# Vd is the dead volume between the valve and the barrel in m^3
	
	return A * x + Vd

def x_2(Pb, Tb, x_dot, Pf, k, MW, Rbar, Patm, Tatm, A, m):
	# returns the current projectile acceleration in m/s^2
	# Pb is the current barrel pressure in Pa
	# Tb is the current barrel temperature in K
	# x_dot is last iteration's velocity in m/s
	# Pf is the minimum pressure to move the projectile in Pa
	# k is the ratio of specific heats of the gas
	# MW is the molecular weight of the gas in kg/kmol
	# Rbar is the gas constant in J/kmol*K
	# Patm is the atmospheric pressure in Pa
	# Tatm is the atmopsheric temperature in K
	# A is the cross-sectional area of the barrel in m^2
	# m is the projectile mass in kg
	# Uses basic Newtonian mechanics (i.e. F = m * a)
	
	# find the total pressures (sum of the static and dynamics pressures)
	
	Pb_t = Pt(Pb, k, x_dot, MW, Tb, Rbar) # barrel gas
	Pa_t = Pt(Patm, k, x_dot, MW, Tatm, Rbar) # atmospheric gas
	
	# calculate the total pressure and friction forces
	
	force = A * (Pb_t - Pa_t)
	friction = A * Pf
	
	# use some logic to determine if the projectile will accelerate and if
	# so, how the force of pressure and friction combine
	
	if x_dot == 0:
		# if the projectile is static and if the total friction force
		# is greater than the total pressure force, then the projectile
		# will not accelerate because the static friction force has not
		# been exceeded
		
		if abs(force) > abs(friction):
			return (force - friction) / m
		else:
			return 0
	else:
		# if the velocity is negative, the direction of the friction
		# force is opposite of the normal direction
		
		if x_dot >= 0:
			return (force - friction) / m
		else:
			return (force + friction) / m

def Tt(m_dot_in, m_cv, T_flow, T_cv, P_cv, work, h, Cv):
	# returns the current barrel or chamber gas temperature in K depending
	# on the direction of flow neglecting the kinetic energy of the flow
	# m_dot_in is the mass flow rate into the control volume in kg/s
	# m_cv is the total mass of the control volume gas at time t - deltat
	# in kg/s
	# T_flow is the temperature of the flowing gas in K
	# T_cv is the temperature of the control volume gas at time t - deltat
	# in K
	# P_cv is the pressure of the control volume gas at time t - deltat in
	# Pa
	# work is the work done during the last time interval in J
	# h is the step size in s
	# Cv is the specific heat capacity at constant volume of the gas
	# Derived from the Reynolds transport theorem
	# See http://en.wikipedia.org/wiki/Reynolds_transport_theorem
	
	return (m_dot_in * T_flow * h + m_cv * T_cv) / (m_dot_in * h + m_cv) - work / (Cv * (m_dot_in * h + m_cv))

def PE(Pstart, Pmin, Vstart, k):
	# returns the potential energy of a gas in J
	# Pstart is the maximum (starting) pressure of the gas in Pa
	# Pmin is the pressure the gas will be expanded to in Pa
	# Vstart is the starting volume of the gas in m^3
	# k is the ratio of specific heats of the gas
	# Derived from polytropic process relationships
	# See http://en.wikipedia.org/wiki/Polytropic_process
	
	return (Pmin * pow(Pstart / Pmin, 1 / k) * Vstart - Pstart * Vstart) / (1 - k)

def Ps(m, T, V, Rbar, MW):
	# returns the pressure of a gas according to the ideal gas law and the
	# conditions given
	# m is the gas mass in kg
	# T is the gas temperature in K
	# V is the volume in m^3
	# Rbar is the gas constant in J/kmol*K
	# MW is the molecular weight of the gas in kg/kmol
	# Derived from the mass form of the ideal gas law
	# See http://en.wikipedia.org/wiki/Ideal_gas_law
	
	return (m * Rbar * T) / (V * MW)

def loop():
	# this is the main loop which is run until the projectile leaves the
	# barrel
	global i, Pb, Pc, Tb, Tc, k, t, to, Cva, SG, MW, h, Cv, Cp, x, A, Vd, Rbar, mb, mc, x_dot_1, x_dot, Tatm, Patm, d, L, m, P, Pf, to, Vc, Vd, Vp, Pr, Rbar, g, Tr, pr, work
	
	if Pb <= Pc:
		# flow into the barrel
		m_dot_in = m_dot(Pc, Pb, Tc, k, t, to, Cva, SG, MW, pr)
		Tb = Tt(m_dot_in, mb, Tc, Tb, Pb, work, h, Cv)
		mb = mb + m_dot_in * h
		mc = mc - m_dot_in * h
	else:
		# flow into the gas chamber
		m_dot_in = m_dot(Pb, Pc, Tb, k, t, to, Cva, SG, MW, pr)
		Tc = Tt(m_dot_in, mc, Tb, Tc, Pc, 0, h, Cv)
		Tb = Tt(-m_dot_in, mb, Tc, Tb, Pb, work, h, Cv)
		mb = mb - m_dot_in * h
		mc = mc + m_dot_in * h
	
	# calculate the new chamber and barrel pressures
	
	Pb = Ps(mb, Tb, Vb(x, A, Vd), Rbar, MW)
	Pc = Ps(mc, Tc, Vc, Rbar, MW)
	
	# save the current velocity for the work calculation
	
	x_dot_1 = x_dot
	
	# Computer the acceleration, then use the forward Euler method to
	# computer the new velocity
	
	accel = x_2(Pb, Tb, x_dot, Pf, k, MW, Rbar, Patm, Tatm, A, m)
	x_dot = x_dot + h * accel
	
	# the increase in kinetic energy between the last step and now is the
	# work done on the projectile
	
	work = m * (pow(x_dot, 2) - pow(x_dot_1, 2)) / 2
	
	# move the projectile down the barrel
	
	x = x + h * x_dot
	
	# increment the time counter
	
	t = t + h

# get gas properties

Cp = gasprop('Cp', gas)
Cv = gasprop('Cv', gas)
k = round(Cp / Cv, 3)
MW = gasprop('MW', gas)
MWair = gasprop('MW', 'air')

# calculate some shortcuts

SG = MW / MWair # calculate the specific gravity of the gas from known values
A = (pi / 4) * pow(d, 2) # calculate the cross-sectional area of the barrel
pr = prcrit(k) # calculate the critical pressure ratio

# initial conditions before looping

t = 0
x = 0
work = 0
x_dot = 0
Pb = Patm
Pc = P + Patm
mb = (Pb * Vb(x, A, Vd) * MW) / (Rbar * Tb)
mc = (Pc * Vc * MW) / (Rbar * Tc)
i = 0

# run the loop

while x < L and t < 0.1:
	loop()
	if x > 0.015 and i == 0:
		print t, x
		i = 1

print "Vf =", round(x_dot, 1), "m/s"
print "tf =", round(t * 1e3, 2), "ms"
print "Tb =", round(Tb - 273.15, 1), "C"
print "a =", round(a(k, MW, Tb, Rbar), 1), "m/s"
print "M =", round(x_dot / a(k, MW, Tb, Rbar), 3)
#print "mair =", 1e5 * (P * Vc * MW) / (Rbar * Tatm)