Decompression Scenarios

from math import *


#Calculation script for decompression scenarios.  See:
#   Description: https://plus.google.com/+IanMallett/posts/246QxPASNrU
#   Calculation addendum: https://plus.google.com/+IanMallett/posts/Jawkrf6C2mM


#Physics parameters
#"μ" is the gas's molar mass: 0.02897 kg mol⁻¹ (unused)
#"R" is the gas constant: 8.3144598 J mol⁻¹ K⁻¹ (unused)
v_a = 343.15

#Human parameters
A = 0.7
C_D = 1.2
m = 65.57

#Chamber parameters
diameter_chamber = 2.0
diameter_hole = 0.1 #2.0
h = pi * (diameter_hole*0.5)**2.0
H = h / (pi * (diameter_chamber*0.5)**2.0)
rho_0 = 1.225
T = 273.15 + 20.0
V = 30.0

#Want to calculate (as applicable)
#f(t)
#v(t)
#a(t)
#M(t)
#rho(t)
#delta-v



#Scenario 1:
def f(t):
    c = (rho_0 * C_D * A) / (2.0 * m)
    return m * c / (( c*t + 1.0/v_a )**2.0)
def v(t):
    c = (rho_0 * C_D * A) / (2.0 * m)
    return (c * v_a * t) / (c*t + 1.0/v_a)
def a(t):
    c = (rho_0 * C_D * A) / (2.0 * m)
    result = c / (( c*t + 1.0/v_a )**2.0)
    return result, result/9.80665
##print(f(0))
##print(v(0))
##print(a(0))
##print()
##print(f(1))
##print(v(1))
##print(a(1))



#Scenario 2 (be sure to set "h" correctly):
def f(t):
    c = (rho_0 * C_D * A) / (2.0 * m)
    return m * c / (( c*t + 1.0/(v_a*H) )**2.0)
def v(t):
    c = (rho_0 * C_D * A) / (2.0 * m)
    return (c * v_a*H * t) / (c*t + 1.0/(v_a*H))
def a(t):
    c = (rho_0 * C_D * A) / (2.0 * m)
    result = c / (( c*t + 1.0/(v_a*H) )**2.0)
    return result, result/9.80665
##print(f(0))
##print(v(0))
##print(a(0))
##print()
##print(f(1))
##print(v(1))
##print(a(1))



#Scenario 3 (be sure to set "h" correctly):
def f(t):
    a = (rho_0 * C_D * A) / (2.0 * m)
    b = -h * v_a / V
    return m * (a * exp(b*t)) / (( (a/b)*(exp(b*t)-1.0) + 1.0/v_a )**2.0)
def v(t):
    a = (rho_0 * C_D * A) / (2.0 * m)
    b = -h * v_a / V
    return v_a - ( (a/b)*(exp(b*t)-1.0) + 1.0/v_a )**-1.0
def a(t):
    a = (rho_0 * C_D * A) / (2.0 * m)
    b = -h * v_a / V
    result = (a * exp(b*t)) / (( (a/b)*(exp(b*t)-1.0) + 1.0/v_a )**2.0)
    return result, result/9.80665
def M(t):
    b = -h * v_a / V
    return V * rho_0 * exp(b * t)
def rho(t):
    b = -h * v_a / V
    return rho_0 * exp(b * t)
def delta_v():
    a = (rho_0 * C_D * A) / (2.0 * m)
    b = -h * v_a / V
    return v_a - ( 1.0/v_a - a/b )**-1.0
##print(f(0))
##print(v(0))
##print(a(0))
##print(M(0))
##print(rho(0))
##print()
##print(f(1))
##print(v(1))
##print(a(1))
##print(M(1))
##print(rho(1))
##print()
##print(delta_v())
##t=0.0; step=0.001 #ad-hoc binary search to find time for various ratios (0.99 below)
##r = v(0.0)/delta_v()
##while step>1e-9:
##    t += step
##    r = v(t)/delta_v()
##    if r >= 0.99:
##        t -= step
##        step *= 0.5
##print((t,r))



#Scenario 4 (be sure to set "h" correctly):
def f(t):
    a = (rho_0 * C_D * A) / (2.0 * m)
    b = -h * v_a / V
    return m * (a * exp(b*t)) / (( (a/b)*(exp(b*t)-1.0) + 1.0/(H*v_a) )**2.0)
def v(t):
    a = (rho_0 * C_D * A) / (2.0 * m)
    b = -h * v_a / V
    return (H*v_a) - ( (a/b)*(exp(b*t)-1.0) + 1.0/(H*v_a) )**-1.0
def a(t):
    a = (rho_0 * C_D * A) / (2.0 * m)
    b = -h * v_a / V
    result = (a * exp(b*t)) / (( (a/b)*(exp(b*t)-1.0) + 1.0/(H*v_a) )**2.0)
    return result, result/9.80665
def M(t):
    b = -h * v_a / V
    return V * rho_0 * exp(b * t)
def rho(t):
    b = -h * v_a / V
    return rho_0 * exp(b * t)
def delta_v():
    a = (rho_0 * C_D * A) / (2.0 * m)
    b = -h * v_a / V
    return (H*v_a) - ( 1.0/(H*v_a) - a/b )**-1.0
##print(f(0))
##print(v(0))
##print(a(0))
##print(M(0))
##print(rho(0))
##print()
##print(f(1))
##print(v(1))
##print(a(1))
##print(M(1))
##print(rho(1))
##print()
##print(delta_v())