# BAGS - http://trettel.org/bags/
# btrettel's air gun simulation
# bags.py - standalone script
# Last updated 05/01/09
# version 1.0.0

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 # interval in seconds

# atmospheric conditions and gun gas temperatures
# available from http://weather.gov in the USA
# 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 = Tr # 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
	
	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
	
	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
	
	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
	
	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
	
	# 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
	
	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
	
	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
	
	return (m * Rbar * T) / (V * MW)

def loop():
	# this is the main loop which is run until the projectile leaves the barrel
	global 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)

# run the loop

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

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)