{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "######The latest version of this IPython notebook is available at [https://github.com/jckantor/Airbag-Design-for-Cargo-Airdrop](https://github.com/jckantor/Airbag-Design-for-Cargo-Airdrop) under the [MIT License](https://github.com/jckantor/Airbag-Design-for-Cargo-Airdrop/blob/master/LICENSE).\n",
    "\n",
    "J.C. Kantor (Kantor.1@nd.edu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Venting Mass Flowrate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mass Flow from an Orifice\n",
    "\n",
    "The mass flow through the airbag vent orifice is given by\n",
    "\n",
    "$$\\dot{m} = C_DA_o \\rho_o v_o$$\n",
    "\n",
    "where the subscript $o$ refers to conditions at the orifice. We begin by constructing an estimate for $\\rho_o v_o$ using Bernoulli's equation for compressible flow on a streamline, the pressure-density relationship for the isentropic expansion of an ideal gas, and an assumption regarding pressure at the orifice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,\n",
       "    900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Crimson+Text' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Kameron' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Lato:200' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Lato:300' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Lato:400' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro' rel='stylesheet' type='text/css'>\n",
       "\n",
       "<style>\n",
       "\n",
       "@font-face {\n",
       "    font-family: \"Computer Modern\";\n",
       "    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "}\n",
       "\n",
       "#notebook_panel { /* establish main background and color */\n",
       "    background: #fff;\n",
       "    color: #ff0000;\n",
       "}\n",
       "\n",
       "div.cell { /* set cell width */\n",
       "    width:800px;\n",
       "    margin-left:0%;\n",
       "    margin-right:auto;\n",
       "}\n",
       "\n",
       "div #notebook { /* centers the content */\n",
       "    background: #fff;\n",
       "    width: 1000px;\n",
       "    margin: auto;\n",
       "    padding-left: 1em;\n",
       "}\n",
       "\n",
       "#notebook li { /* more space between bullet points */\n",
       "    margin-top:0.8em;\n",
       "}\n",
       "\n",
       "/* draw border around running cells */\n",
       "div.cell.border-box-sizing.code_cell.running {\n",
       "    border: 3px solid #111;\n",
       "}\n",
       "\n",
       "/* Put a solid color box around each code cell and its output,\n",
       "     visually linking them together */\n",
       "div.cell.code_cell {\n",
       "    background-color: /*#fff*/ rgba(171,165,131,0.2);\n",
       "    border-radius: 10px; /* rounded borders */\n",
       "    padding: 1em;\n",
       "    margin-top: 1em;\n",
       "}\n",
       "\n",
       ".CodeMirror {\n",
       "    font-family: \"Source Code Pro\";\n",
       "    font-size: 90%;\n",
       "}\n",
       "\n",
       "div.text_cell_render code {\n",
       "    font-family: \"Source Code Pro\";\n",
       "    font-size: 90%;\n",
       "    color: red;\n",
       "}\n",
       "\n",
       "pre {\n",
       "    font-family: \"Source Code Pro\";\n",
       "    font-size: 90%;\n",
       "}\n",
       "\n",
       "div.text_cell_render { /* rendering of basic text cell */\n",
       "    font-family: 'Alegreya Sans', Computer Modern, \"Helvetica Neue\", Arial, Helvetica, sans-serif;\n",
       "    line-height: 135%;\n",
       "    font-size: 120%;\n",
       "    width:700px;\n",
       "    margin-left:auto;\n",
       "    margin-right:auto;\n",
       "}\n",
       "\n",
       ".text_cell_render h1 {\n",
       "    font-family: 'Alegreya Sans', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 32pt;\n",
       "    line-height: 100%;\n",
       "    color:#CD2305;\n",
       "    margin-bottom: 0.4em;\n",
       "    margin-top: 0em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       "h2 {\n",
       "    font-family: 'Fenix', serif;\n",
       "}\n",
       "\n",
       "h3 {\n",
       "    font-family: 'Fenix', serif;\n",
       "    margin-top:12px;\n",
       "    margin-bottom: 3px;\n",
       "}\n",
       "\n",
       "h4 {\n",
       "    font-family: 'Fenix', serif;\n",
       "}\n",
       "\n",
       ".text_cell_render h5 {\n",
       "    font-family: 'Alegreya Sans', sans-serif;\n",
       "    font-weight: 100;\n",
       "    font-size: 16pt;\n",
       "    color: #CD2305;\n",
       "    font-style: italic;\n",
       "    margin-bottom: .5em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h6 { /*use this for copyright note*/\n",
       "    font-family: 'Alegreya Sans', sans-serif;\n",
       "    font-weight: 400;\n",
       "    font-size: 10pt;\n",
       "    font-style: normal;\n",
       "    line-height: 120%;\n",
       "    margin-bottom: 1px;\n",
       "    margin-top: 1px;\n",
       "}\n",
       "\n",
       ".MathJax, .MathJax_Display  {\n",
       "    font-family: 'Alegreya Sans', sans-serif  !important;\n",
       "    font-size: 90%  !important;\n",
       "    vertical-align: 0px  !important;\n",
       "}\n",
       "\n",
       ".warning{\n",
       "    color: rgb( 240, 20, 20 )\n",
       "}\n",
       "\n",
       "</style>\n",
       "\n",
       "<script>\n",
       "MathJax.Hub.Config({\n",
       "    TeX: {\n",
       "        extensions: [\"AMSmath.js\"]\n",
       "    },\n",
       "    tex2jax: {\n",
       "        inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "        displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ],\n",
       "        processEscapes: true\n",
       "    },\n",
       "    displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "    \"HTML-CSS\": {\n",
       "        styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "    }\n",
       "});\n",
       "</script>\n",
       "\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Initializations\n",
    "from IPython.core.display import HTML\n",
    "HTML(open(\"../styles/custom.css\", \"r\").read())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import math\n",
    "pi = math.pi\n",
    "\n",
    "from pint import UnitRegistry\n",
    "ur = UnitRegistry()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Coefficient of Discharge\n",
    "\n",
    "Based on Browning (1963)\n",
    "\n",
    "$$C_D = 0.9 - 0.3\\frac{P_a}{P}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Flow under Unchoked Conditions\n",
    "\n",
    "### Bernoulli's Equation\n",
    "\n",
    "Bernoulli's equation for the behavior of a compressible fluid along a streamline\n",
    "\n",
    "$$\\frac{v_o^2}{2} + \\frac{\\gamma}{\\gamma-1} \\frac{P_o}{\\rho_o}  = \\frac{\\gamma}{\\gamma-1} \\frac{P}{\\rho}$$\n",
    "\n",
    "which, for a streamline starting deep within the airbag, gives for $v_o$ \n",
    "\n",
    "$$v_o = \\sqrt{\\frac{2\\gamma}{\\gamma -1}\\left(\\frac{P}{\\rho} - \\frac{P_o}{\\rho_o}\\right)}$$\n",
    "\n",
    "then for $\\dot{m} = C_D A_o\\rho_ov_o$\n",
    "\n",
    "$$\\dot{m} = C_D A_o \\sqrt{\\frac{2\\gamma\\rho_o^2}{\\gamma -1}\\left(\\frac{P}{\\rho} - \\frac{P_o}{\\rho_o}\\right)}$$\n",
    "\n",
    "### Adiabatic Expansion\n",
    "\n",
    "For the adiabatic expansion of an ideal gas,\n",
    "\n",
    "$$\\frac{\\rho_o}{\\rho} = \\left(\\frac{P_o}{P}\\right)^{\\frac{1}{\\gamma}}$$\n",
    "\n",
    "Substituting into the expression for $\\rho_ov_o$,\n",
    "\n",
    "$$\\dot{m} = C_D A_o \\sqrt{\\frac{2\\gamma\\rho^2}{\\gamma -1}\\left(\\frac{P_o}{P}\\right)^\\frac{2}{\\gamma}\\left(\\frac{P}{\\rho} - \\frac{P_o}{\\rho\\left(\\frac{P_o}{P}\\right)^\\frac{1}{\\gamma}}\\right)}$$\n",
    "\n",
    "After rearrangement\n",
    "\n",
    "$$\\dot{m} = C_D A_o \\sqrt{\\frac{2\\gamma\\rho P}{\\gamma -1}\\left[\\left(\\frac{P_o}{P}\\right)^\\frac{2}{\\gamma} - \\left(\\frac{P_o}{P}\\right)^{1+\\frac{1}{\\gamma}}\\right]}$$\n",
    "\n",
    "### Evaluating at the Orifice\n",
    "\n",
    "Assuming the pressure at the orifice is equal to the surrounding ambient pressure, $P_o = P_a$.  With that assumption we have\n",
    "\n",
    "$$\\dot{m} = C_D A_o \\sqrt{\\frac{2\\gamma\\rho P}{\\gamma -1}\\left[\\left(\\frac{P_a}{P}\\right)^\\frac{2}{\\gamma} - \\left(\\frac{P_a}{P}\\right)^{1+\\frac{1}{\\gamma}}\\right]}$$\n",
    "\n",
    "This expression is found in [_Handbook of Chemical Hazard Analysis Procedures_](http://nepis.epa.gov/Exe/ZyNET.exe/10003MK5.TXT?ZyActionD=ZyDocument&Client=EPA&Index=1986+Thru+1990&Docs=&Query=&Time=&EndTime=&SearchMethod=1&TocRestrict=n&Toc=&TocEntry=&QField=&QFieldYear=&QFieldMonth=&QFieldDay=&IntQFieldOp=0&ExtQFieldOp=0&XmlQuery=&File=D%3A%5Czyfiles%5CIndex%20Data%5C86thru90%5CTxt%5C00000003%5C10003MK5.txt&User=ANONYMOUS&Password=anonymous&SortMethod=h%7C-&MaximumDocuments=1&FuzzyDegree=0&ImageQuality=r75g8/r75g8/x150y150g16/i425&Display=p%7Cf&DefSeekPage=x&SearchBack=ZyActionL&Back=ZyActionS&BackDesc=Results%20page&MaximumPages=1&ZyEntry=1&SeekPage=x&ZyPURL), Appendix B, Federal Emergency Management Agency, U.S. Dept. of Transportation, and U.S. Environmental Protection Agency, 1989. The relevant section is on page B-4 located on the page 391 of the pdf file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x102b341d0>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEZCAYAAABiu9n+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmUHGW5x/HvLwkhrIZcJCQQVtkiKGtYvQSRyA4qSxCQ\nTdSLAipe2VQiIiJHFhfcQLgBMRhZIsgalGFTFhEQCBgWgwRMIMoOQiDP/eN9J+lMpmd6Jl1d3TO/\nzzl9UlVdXfV0T6efetdSRGBmZtaZAWUHYGZmzctJwszMqnKSMDOzqpwkzMysKicJMzOryknCzMyq\ncpLooyStIWmepB7/jSWNlfRMneKYJ2mtRr2uB8dfStI1kl6S9OuiztPMJB0o6caK9W0lPS7pFUl7\nSbpO0qfKjNHK5yRRss5+DCVNkHRJWTE1iqQ2SW9KerXisWWDTr8PsBIwLCL2b9A5CyNpd0n3SHpN\n0hxJv5S0SleviYhLI+KjFZtOBX4QEctHxG8jYteIuLjYyEHSJyX9Of/9n8vJadv83ARJc3PiekXS\n3yT9UNLKBcSxoaQbJb0gaV69j9+qnCSaU38Z4RjA5yNiuYrH3Q069+rA9Ijo9MdA0sAGxbHYJO0D\nXAqcDfwX8H7gLeAOSUOrvKaz97caMK2oOKvE8WXgHOA0UtIeBZwH7FGx26SIWB5YAfgYsDJwXwGJ\n4m3gMuCIOh+3pTlJNCfNX0hVPzMlfVnS7HyldWjF80tJOkvSjFx1crukJSuOdZCkp/PV0UkVr1tS\n0rmSns2PcyQN7jQY6RhJj0gamV/3vXzMWZJ+ImlIxb7/m2OcKenwunwY0nskXSzp+fw+T5ak/NzT\nkjbNywfmktkGef0ISVd1crxvAl8H9s9Xr4dLOlTSnZLOljQHOEXS8l2ct3L/FyU9IWkbSYdJ+kf+\nW1Wtqsn7TctXx09K+kzFc13+zTscR8BZwLci4rKIeCsiZgOfBl4DvtRJvHOACXnb7fn5J4G1gGty\nTIOVSnpHVJzryIqYH5G0Sd4+UtIV+XN6StLRtf5dgW8CR0XElIh4MyLejYhrI+KEyl0B8nPTgP2B\nF4DjajlPrSJiekRcRIMTZbNzkmgNw4HlgZGkq5zz8n8wgO8BmwBbA8OA/2Xhksi2wLrAjsA3JK2X\nt58MjAE+mB9jgK91PLGkbwCfAv47Ip4DzgDel1/zPmAV4Bt5351J/3E/ks/5kRrem7rfhR8CywFr\nAtvneA7Lz7UBY/Py9sCT+d/29baOB4uIU4DTgcty6eXC/NSY/PqV8vM/6uK87fs/SPrcJwGTgU2B\ntYGDgB9JWrrKe5oN7JavkA8Dzmn/0c26+ptXWo909f2bDu8xgCuAnTrE2/7+vt1h/7WBfwC75+qm\nt0nfowCQtC9wCnBwjnlP4F9KbV7XAPfnWHcEvihpXJX3XWlrYAiwSCLvSi79/Rb4UGfPS9ouJ+5q\nj216cr5+LyL8KPEBzAPW6rBtAnBJXh4LvAEMqHh+Nuk//ID83EadHHeNfOyRFdvuBvbLy08AO1c8\nNw74e8U5Z5KqL24DlsvbRbo6XavidVsDT+XlC4HTK55bp7P3V/F8G/A68GJ+/Lnj5wIMJFWdrF/x\n3GeAW/Ly4cBv8/K0vD4pr88ANq5y7vmfcV4/FHi6Yr278x5Kqq5qf26jHPN7K7bNAT5Q4/fgKuCY\n7v7mnbxuu3zewZ0897n2GDu+v4ptt1es/x34cMX6LcDheflG4OhOzrFlJ8c9Ebiwhvd8IPDPbvZZ\n6O/U2Xur94N08TOviGO34mMQVrZ3gSU6bFsCmFux/q9YuO78DWBZYEXSldiTXRx/Vievg3TV93TF\nc//I29oNJVVZjI+IV/O29wJLk+qD2/cTC0qkI4B7OxyzK0H64bmwi31WJH0eHWNtb5S9DfieUv30\nQNIV9QRJqwPviYgHuomhUmWPru7OC+mHu92bABHxQodty9IJSbuQrszXIX1+SwN/rdil2t+8ozn5\n3xEdYm3fVhnP4vRYW5XOv2erAyMlvVixbSDp79KdfwErShoQVdqGurBKfr0VzNVN5fsHqTqj0pqk\nq+DuzAH+Q7ry6annSKWNdqvlbe1eBHYHLqoons8h/fCNjogV8mNopOoHgH/m41Qec3HNISXMjrHO\nBIiIJ0g/oEcDt+aENot01X97F8ftrHNA5bYuz7s4cpvRFcCZwEoRsQJwHbVVvXX0txzTfh3OMQD4\nBPD7is2L0yHiGTr/nv2DVAJdoeKxfETsXsMx/0QqrX2si30WiTm/tz2o8veV9CEt3GOu42PbGmKz\nzEmifL8GviZpFUkDJH2E9ON8eXcvzFdfFwJnSxohaaCkrVWlAbqDSfm8K0pakdSusFC324i4jVQl\ncKWkLfL5zgfOlfRegBx3e/3zZOBQSRvkuvhTaoijyx/GiHg3H/fbkpbNJYQvAb+s2O1W4Av5X0jV\nWJXrRZ23twbnxxxgXi5V1FKH31mcAXyF9Lc8QNKQXKq6gFTyOKcO8ZKP9xVJmyp5n6TVgHuAVyV9\nVakTxUClrqSbw/xG+E5LCRHxMul7d57SuIylJS0haRdJ3827VXbiGKTUKWESqV3l7CrHvT0W7jHX\n8XFntTep1AljcF5eUgt3AumXnCTKdyrwR+AO4N+khuFPRurF0a6rK8CvAA+Rqnn+BXyHBf+xunrd\nacCfSVUcf83Lp3U8Z0TcTKrnv0bSxsDxpPaMuyS9DEwlNVITETcA5wJ/AKaTrmK7u3qt9nzl9qNJ\nbRdPka4eLwUuqnj+VtIP4m1V1qsdP7pY7+68ne1f05V6Lu0cQ0pC/wYOIDXE9vhY+XiTgYNJSWwO\n8AiwJLBtRLRXA1WLt9aYLyc1dv8KeAW4ElghXzjsDmxM+pxeAH5OanSH1Khe9Uc5Is4GvkzqNPE8\nqWRyFAsas4PcCw14ifQ5vQBsFhGzFj1i70lag1QqfTif903g0XqeoxUpN9QUc/DUk+ayik1rkboe\n/pJ0Bb06qVplv4h4Kb/mRNKP0rukhrybCgvQzAol6XxgckRMLTsW651Ck8RCJ0r1iM+SeuUcDcyJ\niDMlHU+6IjlB0mjSlcoWpIapm4F1e9GoZWZmddDI6qaPAE9ExDOkPtYT8/aJwN55eS9S98W5ETGD\nVK0xpoExmplZhUYmifGkBieA4ZFGhULqRjg8L49k4d4jM1m4y6GZmTVQQ5JE7m2zBx1GhcL83hld\n1Xn1l3mMzMyaTqMG0+0C3Fcx0Gi2pJUjYpakEaReDZDaLEZVvG7VvG0+SU4aZma9EBE9HovTqCRx\nAAuqmgCuBg4Bvpv/nVKx/VeSziZVM61D6oe9kN68UeucpAkRMaHsOPoKf57148+yvnp7gV14kpC0\nDKnR+siKzWcAk/MMkzPIo0UjYpqkyaQ5eN4hzQ7pkoOZWUkKTxIR8TppHpzKbf+mygyhEXE6aQZO\nMzMrmUdcW1vZAfQxbWUH0Ie0lR2ANXAwXb1ICrdJmJn1TG9/O12SMDOzqpwkzMysKicJMzOryknC\nzMyqcpIwM7OqnCTMzKwqJwkzsz5O6v3AaScJM7M+SGJViSMkLifd8rVXnCTMzPoAicESO0icKfEQ\n8CBp+qNrgNG9Pq5HXJuZtSaJ1Ui3YtgF2AH4G3B9ftwbwbsL9u3db6eThJlZi5AYDGwL7EpKDMOB\nG0lJ4aaI6tVKThJmZn2QxCqkhLAr8GFgOnAdKTH8ubK00PVxnCTMzFqexEBgS2A3UmJYHbgJuBa4\nMWL+nTx7eFwnCTOzliQxDPgoKTHsTLpt87WkEsNdEbyz+OdwkjAzawkSIvU42p2UGDYm3T/jWuC6\nCJ6p/zmdJMzMmpbEEGB7UmLYHRDwu/xoi+A/xZ6/d7+dhd++1Mysv5IYTiop7EFqdH6ING5hd2Ba\nBE1/le6ShJlZneRqpA+QksLuwHrAVFJiuD6COeXF5uomM7OGk1iSVI20Jyk5vAtcTUoMt0fwdonh\nzefqJjOzBpFYgdQ9dS9gHDCNlBh2pUWqkWrlkoSZWQ0k1iQlhT2BzUi9kX4L/K63Yxcaqbe/nYVP\n8CdpqKTLJT0qaZqkLSUNkzRV0nRJN0kaWrH/iZIel/SYpHFFx2dm1hkJSWwu8S2JvwJ3ARsC5wAj\nItgrggtbIUEsjsJLEpImArdGxIWSBgHLACcDcyLiTEnHAytExAmSRgO/ArYAVgFuBtaNiHkVx3NJ\nwswKIbEEqX1hb1Kp4Q1gCqnEcHetU2A0o6ZsuJb0HuD+iFirw/bHgO0jYraklYG2iFhf0onAvIj4\nbt7vBmBCRNxV8VonCTOrG4llSaOc9ya1KUwnJYYpETxWZmz11KwN12sCL0i6CPggcB/wRWB4RMzO\n+8wmzWQIMJJUpGs3k1SiMDOrG4kVST2RPgaMBf4EXAV8NYLnSgyt6RSdJAYBmwJfiIh7JZ0LnFC5\nQ0SEpK6KM4s8J2lCxWpbRLTVIVYz68MkRpFKCx8n/S5NBX4NfCqCl8qMrQiSxpIS4GIpOknMBGZG\nxL15/XLgRGCWpJUjYpakETC/4edZYFTF61fN2xYSEROKC9nM+gqJdYBPkBLDWqQpMM4BpkbwZpmx\nFS1fPLe1r0s6pTfHKbR3U0TMAp6RtG7e9BHgEdIgk0PytkNI9X+Q+hmPlzRY0prAOsA9RcZoZn1H\n7pG0kcSEfAvPW0kXnieReiQdGsHVfT1B1FMjBtMdDVwqaTDwJHAYMBCYLOkIYAawH0BETJM0mTQw\n5R3gqGi1gRxm1lB5KozNgH1IpYYlgCuAzwF/imBeFy+3bngwnZm1HIkBwBhSYtgHeJuUGK4A7utL\nI57rpVl7N5mZ1UVODFsD+5JKDK+QksIewMNODMWomiQkbUbqWSQ66WEUEX8pMC4zs/bEsC0LEsO/\ngd8AH41gWpmx9RdVq5sktZGSw1Kk+r6/5qc+APw5IrZuRICdxOXqJrM+rKLEsB+pKulfpMTwm740\nuK3R6l7dFBFj84GvBI6MiIfy+obAN3sZp5nZInLj81akxLAv8CIwGdjRiaFctbRJrN+eIAAi4mFJ\nGxQYk5n1AzkxbA7sT0oMr5MGt+0UwaNlxmYL1JIk/irpAuCXpPaJTwIPFhqVmfVJOTFsBIwnJYd3\nSYlhtwgeLjM261y3XWAlLQX8D/ChvOk24CcRUehNu7uIx20SZi1GYl1SYhhPmgn6MlJyuN+9khqj\nKWeBLYKThFlryHMljQcOAEaQ2hh+TRrg1lo/PH1AYUlC0nbAKcAaLKieio7TfzeKk4RZ88qzq+5L\nqpYeDVwJTAJubeV7MfQFRSaJv5Gm9/4LLPgjR8Scnp6sHpwkzJpLvh/DXqTEsB1wPenmYTdG8FaZ\nsdkCRY64fikiru9FTGbWR+U7uI0DDiTdqOePwKXA/hG8VmZsVl+1lCTOIE3IdyUsuCooa8S1SxJm\n5cg9k7YEDiKNZ3iClBgmR/BCmbFZ94osSWxFGnm9eYftO/T0ZGbWeiTeR0oMB5FmZ/4lsFUET5Ua\nmDWEezeZ2SIk/os0juFg0m2ILyMlB8+w2qKKbLgeSurd9N95UxtwakS83NOT1YOThFkxJJYEdiMl\nhh2A64BLSHdxe6fM2GzxFZkkrgQeAiaSRlwfDHwgIj7em0AXl5OEWf1UtDN8itTO8DDp//oVEbxS\nZmxWX0UmiQcj4oPdbWsUJwmzxZcHuh1Mun2wgIuBSyJ4utTArDBFNly/KelDEXF7PtF2wBs9PZGZ\nlUtiaeBjwKHApqTptw8F7nI7g1VTS5L4HHCxpPfk9RdJVx9m1uRyddI2pGTwCeAu4ALgtxGUMv+a\ntZaaeze1J4myGqwr4nB1k1k3JEaS2hkOI3Vh/z9SddKzZcZl5entb+eAGg78HUlDI+LliHhZ0gqS\nTutdmGZWFInBEp+QuJbUAL0WqQSxQQRnOEFYb9TScP1ARGzcYdv9EbFJoZFVj8clCbMKEhsCR5Cm\nyHgEuBC4MoLXSw3MmkqRDdcDJA1pv39Evr/E4J6eyMzqR2I50mC3TwOjSNVJ20TwRJlxWd9TS5K4\nFPi9pAtJXeUOI3WXq4mkGcArpBlk50bEGEnDSPPKrw7MAPaLiJfy/icCh+f9j4mIm2p+N2Z9WMV9\noD8NfJw0sPU04AYPdrOi1NRwLWkXYMe8OjUibqz5BNLfgc0i4t8V284E5kTEmZKOB1aIiBMkjSZN\nMbwFsApwM7BuRMyreK2rm6xfkRhGGtNwJLAk8Avg/yKYVWpg1lKKrG4CeBR4JyKmSlpa0nIR8WpP\n4uuwviewfV6eSLoiOoE0J/2kiJgLzJD0BDCG1G3PrN/IpYbtgM8Ae5CmyDgaaPOYBmukWno3fYY0\n6OanedOqwJQenCOAmyX9WdKRedvwiJidl2cDw/PySGBmxWtnkkoUZv2CxDCJL5IaoH8O3AesHcEn\nI7jFCcIarZaSxOepuJqPiOmSVurBObaNiH9Kei8wVdJjlU9GREjq6ou/yHOSJlSstkVEWw/iMWsq\nFQPePksqZf8uL9/hpGC9JWksMHZxj1NLkngrIt6S1H7iQXTyw11NRPwz//uCpKtICWe2pJUjYpak\nEcDzefdnST012q2at3U85oRaz2/WrCSWJ92j4XPAEFJp/UsR/KvUwKxPyBfPbe3rkk7pzXG6rW4C\nbpV0MrC0pJ1IVU/X1HLw9vaLvLwM6XaHDwFXs2Bqj0NYUH11NTBe0mBJawLrAPfU+mbMWoHExhI/\nA54mXel9EVgvgrOdIKzZ1DKYbiBpoM64vOlG4IKooVtU/qG/Kq8OAi6NiO/kLrCTgdVYtAvsSaQu\nsO8Ax3bsSeXeTdaK8r0a9gWOIn3vfwZcEME/Sw3M+o3CpgrvcJJhwKiIeLCnJ6oXJwlrJRKrk6qT\nDgceBH4CXONxDdZoRc7ddKuk5XOCuA84X9I5vQnSrD+QkMRHJKYAfwGWAv47gnERXOUEYa2klobr\n90TEK5I+DVwcEadIeqjowMxaTZ4q4xBSj8C5wI+AAz2HkrWyWhquB+YeSPsB1+Zt7pZnlkmsI/F9\nFjREfw74YAQ/d4KwVldLkjiV1Fj9ZETcI2lt4PFiwzJrbrlKaVyelvtO4HVSYtgngls9vsH6ih41\nXDcDN1xbmSSWIc2jdCzwNvB9YFIEb5YamFk3ip67yaxfkxgFfIHUS+lOUldWz6NkfV4t1U1m/ZbE\nGIlJpO6rSwJbRbC351Gy/qJqkpC0jdrn4jDrRyQGSnxc4g7SoM97gTUj+GIET5YcnllDVW2TkPRT\nYEtgOnA9cENElD5/vdskrCi5veEw0jQZc4CzSbcB9bgGa3mFjbiWtAGwC2lajqHAH4AbgDsj4t1e\nxLpYnCSs3iRWJrU3fBa4DTgrgj+WG5VZfTVqWo6lgR1ISWPriNispydcXE4SVi8S6wNfAT5BuiPi\nOb5HtPVVDUkSzcBJwhaXxLbAV0n3i/4xcF4Ec8qNyqxY7gJr1gWJAcBuwPGkOyB+DzgggjdKDcys\nyTlJWJ8msQQwnpQc5gJnAFe4MdqsNk4S1idJLEW6D8pXgL8DxwE3eWyDWc90NU5iNUmXSbpD0kmS\nlqh4bkq115mVSWJ5iRNIiWEnYHwEO0RwoxOEWc91NeL6QtL9UY8m1eHeKmnF/NzqBcdl1iMSK0p8\nC3gK2Aj4SAR7RXBXyaGZtbSuqpveGxE/zctfkHQQcJukPRoQl1lN8hiH40hVS1eQps1wN1azOukq\nSQySNCQi/gMQEb+UNIs0bfgyDYnOrAqJkaTG6IOBS4GNI/hHuVGZ9T1dVTf9gtSPfL6IuJl0M/eH\niwzKrBqJVSV+RPoOvgO8P4KjnSDMiuHBdNYS8lTdJ5C6s/4C+F4Ez5cblVnr6O1vp6cKt6aWSw7n\nAQ8ArwHrR/BVJwizxnCSsKYkMULiB6T7OLxOSg7HR/BCyaGZ9SuFJwlJAyXdL+mavD5M0lRJ0yXd\nJGloxb4nSnpc0mOSxhUdmzUfiZUkziK1OcwFRueSg5ODWQm6TBKSfqrkx4txjmOBaTB/INMJwNSI\nWBf4fV5H0mhgf2A0sDPwY0ku6fQTEitIfBt4FBgMbBTBcRHMLjk0s36tqxHXqwN3AL8F/pjXe0TS\nqsCuwAVAe4PJnsDEvDwR2Dsv7wVMioi5ETEDeAIY09NzWmuRWFbiZOBxYCVg09xb6bmSQzMzui5J\njAXWII1eXSOv99Q5wP8C8yq2DY+I9qvD2cDwvDwSmFmx30xglV6c01qAxJISx5IuBjYEtongyAie\nLjk0M6tQdTBdREyUdD5prMS3IuK0nhxY0u7A8xFxv6SxVc4Rkrrqg9vpc5ImVKy2RURbT2Kz8kgM\nBA4ETiW1O+wcwQPlRmXW9+Tf3bGLe5zuZoE9JSJmSzqlF8feBthT0q7AEGB5SZcAsyWtHBGzJI2A\n+V0ZnwVGVbx+1bxtERExoRfxWIkkBOwOnA68DBwUwR3lRmXWd+WL57b29V7+jnfbu+kbkgR8vacH\njoiTImJURKxJGgD1h4g4GLgaOCTvdgjQPqPs1cB4SYMlrQmsA9zT0/Na85HYCrgV+A5wEvAhJwiz\n1lBow3UH7VVHZwA7SZoOfDivExHTgMmknlDXA0dFqw0Ht4VIrCNxOfAb4P+AD0ZwjafsNmsdVafl\nkHQIqfrnCNI0CM9ExMROd24gT8vR/CRWBL4BfBI4C/i+bxNqVq66T8uRE8LqpIbr1ZohQVhzkxgi\ncTzwGKnL8wYRfMcJwqx1dTnBn6SREfGcpBER8c8GxlWVSxLNJzdK70eqOnwAOD6C6eVGZWaVevvb\n6VlgbbFIjAHOJfVg+3LEgt4UZtY8PAusNZTEKhIXA1cB5wObO0GY9T3dzd00QNI2jQrGml9udziZ\nNDvrM6TZWS+KWGhUvZn1EV0miYiYByzO5H7WR0hIYm9SF+VNgS0iODmCV0sOzcwK1N2Ia4CbJe0D\nXOFxC/2TxPrAD0jzax0Zwe9LDsnMGqTbhmtJrwFLA+8C/8mbIyKWLzi2avG44bpBJJYjjbY/HDgN\nOC+CueVGZWa90dvfzm5LEhGxbO9CslZV0aX1LNI9PzaMYFa5UZlZGbpNEvnGPwcCa0bEqZJWA1aO\nCM+r1AdJrAucR5rCfbznWDLr32rpAvtjYGvSFAuQbkbvxuw+RmIpiVOBP5LmztrMCcLMamm43jIi\nNpF0P0BE/FvSEgXHZQ0ksRPwE9Jo6Y0jFrr5k5n1Y7UkibclDWxfkfRecJ/4vkBiJdLdA7cFPh/B\ntSWHZGZNppbqph+SRtWuJOl04E7SfQGsReUxD4cBDwHPAe93gjCzztQ0d5OkDYAd8+rvI+LRQqPq\nOhZ3gV0MEu8DfgYMBT4dwf0lh2RmDVDY3E2SLomIRyPiR/nxaL4NqbUQiYESxwF3AdcBWzpBmFl3\nammT2LByRdIgYLNiwrEiSGwIXEjqmbZlBE+WHJKZtYiubl96kqRXgY0kvdr+AJ4n3Y/ampzEEnky\nvluAC4AdnSDMrCdqmZbjjIg4oUHxdMttErWR+ABwEfACab6lZ0oOycxKVOT9JE6SdLCkb+QTrSZp\nTI8jtIaQGCRxImk6jR8DuzhBmFlv1VKS+ClpXMSHI2J9ScOAmyJi80YE2Ek8LklUkafUmAi8ARwW\nwT9KDsnMmkSRJYktI+Io4E1II64Bj7huInncw+dJU2pcCuzkBGFm9eAR1y1OYiSp59IKwDYRTC85\nJDPrQwobcS1piKS7JT0gaZqk7+TtwyRNlTRd0k2Shla85kRJj0t6TNK4Xr6nfkPi48BfSGMftnOC\nMLN6K3TEtaSlI+KNPLbiDuArwJ7AnIg4U9LxwAoRcYKk0cCvgC2AVYCbgXXzLVQrj9nv2yQkliHN\nufRh4KAI7io5JDNrckW2SQDMAm4H/gQsJWnTWl4UEW/kxcHAQOBFUpKYmLdPBPbOy3sBkyJibkTM\nAJ4A3IuqA4lNgPuAJYFNnCDMrEi13HToW8ChwFMs3BaxQw2vHUCqDlkb+ElEPCJpeETMzrvMJt3c\nBtL9kyt/8GaSShTG/LvFHQN8DTg2gl+VHJKZ9QO1NFzvD6wdEW/39OC5qmhjSe8BbpS0Q4fnQ1JX\n9V2dPidpQsVqW0S09TS2ViIxjDQwbiRpWo2nSg7JzJqcpLHA2MU9Ti1J4hFSz5nZ3e1YTUS8LOla\n0pxPsyWtHBGzJI0gTfMB8CwwquJlq+ZtnR1vQm9jaTUS2wCTgMuBfSPocbI2s/4nXzy3ta9LOqU3\nx6mlTeJ04P7cE+ma/Oh27iZJK7b3XJK0FLATcD9p3qdD8m6HAFPy8tXAeEmDJa0JrAP02/to57EP\nx5F6ln0+guOcIMys0WopSVwMnAE8zII2ie67RMEIYGJulxgAXBIRv8+3QZ0s6QhgBrAfQERMkzQZ\nmAa8AxwVtXS96oMkhpKql1YBxkTwdMkhmVk/Vcu0HPdGxBYNiqdbfb0LbJ6Y70rgeuArEbxVckhm\n1gf09rezliRxNvAWqTpo/g9WRPylpyerh76cJCQOBM7FvZfMrM6KTBJtdFK9FBHddoEtQl9MEhJL\nAGcBuwIfi+ChkkMysz6msCTRbPpakpB4L6nn0quk0dMvlRySmfVBRY+4tgJIbAzcS5qyZC8nCDNr\nNrX0brICSOwD/AT4QgS/LjseM7POVE0SkkZGxHONDKY/yNNrnAx8BhgXwf0lh2RmVlVXJYlf5LvQ\n3QLcANwREe80Jqy+SWIIcAGwLml6jX+WHJKZWZe6bLjOI6XHAjsD2wLPkPrv3xARpdz5rFUbrnMD\n9RTSxIWHRqQ7/ZmZNUJDejdJWgvYBfgosHJENHwq71ZMEhLrANcBk4GvR/jOfmbWWA3vAitpyYho\n+GjgVksSEluT5l/6RgQ/LzseM+ufevvb2eveTWUkiFYjsTfwc+CQCK4vOx4zs55yF9iCSHwaOBXY\nJYL7yo7HzKw3nCTqLHdxPQk4HNg+gsdLDsnMrNe6HXEtaQ9J90t6UdKr+fFKI4JrNRIDgHOAfYHt\nnCDMrNXVMsHfk8DHgIfz7UhL1awN1xIDgfOB9YDdPMWGmTWTIhuuZwKPNEOCaFZ5FtdLgBVJo6hf\nLzkkM7NCIgPrAAAOTUlEQVS6qCVJHA9cL+kWmH/7zIiIs4sLq3XkUdSTAQG7R/CfkkMyM6ubWpLE\nt0jTWA8BBhcbTmvJCeJK4HXgkxHMLTkkM7O6qiVJjIiInQqPpMVUJIjXSAnC81qZWZ9Ty/0krpP0\n0cIjaSEVCeJVnCDMrA+rpXfTa8DSpPaI9uqUiIjlC46tWjyl9m6SGMyCKqYDnSDMrBX49qUNOTeD\ngF+R2mb2dRuEmbWKQudukrQCsA6p8RqAiLitpydrZXmg3PnAUGAPJwgz6w9qGXF9JHAbcBPwTeBG\nYEItB5c0StItkh6R9LCkY/L2YZKmSpou6SZJQytec6KkxyU9Jmlcb95UveWpNr4PvA/4WASe3NDM\n+oVaGq6PBcYAMyJiB2AT4OUajz8X+FJEvB/YCvi8pA2AE4CpEbEu8Pu8jqTRwP7AaNKNjn4sqZYY\ni3YysB1pHIQHyplZv1HLD/B/IuJNAElDIuIx0tQT3YqIWRHxQF5+DXgUWAXYE5iYd5sI7J2X9wIm\nRcTciJgBPEFKUKXJs7keTprNtdbkaGbWJ9Q0LUduk5gCTJX0IjCjpyeStAapFHI3MDwiZuenZgPD\n8/JI4K7Kc5OSSikk9iANJtw+glllxWFmVpZuk0REtF/lT5DUBiwP3NCTk0haFrgCODYiXpUWNLBH\nREjqqovVIs9JmlCx2hYRbT2JpxYSWwIXkqqYptf7+GZmRZI0Fhi7uMfpMklIGkSa/XV9gN78GEta\ngpQgLomIKXnzbEkrR8QsSSOA5/P2Z4FRFS9fNW9bSERM6GkcPSGxGmksxBER3F3kuczMipB/r9va\n1yWd0pvjdNkmERHvAH+TtHpvDq5UZPgFMC0izq146mrgkLx8CKkqq337eEmDJa1J6nZ7T2/O3VsS\ny+Y4zo7g6kae28ys2dQy4vp2UlvCPTC/Z09ExJ7dHlzajtR99q8sqDY6MR9rMrAaqX1jv4h4Kb+m\n/a5u75Cqp27scMzCBtPle0JcCbwAHBmxaFWXmVkrKmzEda7XWkQR7QC1KDhJnEHqqjsuYv606GZm\nLc/Tciz2cfkE8D1g8wj+Ve/jm5mVqbBpOSS92snml4F7geMi4qmenrTZSKwP/BTY2QnCzGyBWsZJ\nfB94BpiU18cDawP3k7qIji0ksgaRWI7UDnFCBPeVHY+ZWTOppU3irxHxgQ7bHoiIjSU9GBEfLDTC\nReOpW3VTnpPp18DLERxZj2OamTWj3v521jItxxuS9pc0ID/2g/n3cW6tBo1FfZY0ad/RZQdiZtaM\nailJrE2qctoqb7oL+CJpkNtmEXFHoREuGk9dShISo0ndc7eN4G+LH5mZWfNy76YeHYMlSXNInRfB\n+fWJzMyseRVW3ZTvCXGVpBfy4wpJq/YuzKbxHeBJ4IKyAzEza2a1tElcRJqmYmR+XJO3tSSJccC+\neES1mVm3ammTWKQHUxm9mirO3evqptzd9WHg0xFMrW9kZmbNq8jeTf+SdLCkgZIGSToImNPzEJvC\n6cAfnCDMzGpTS0liDeCHLOjd9Efg6Ij4R6GRVY+ndy30YhvgcmDDCP5d/8jMzJqXezd1+RqGkEaI\nfz2Cy4uJzMysedV97iZJP+zidRERx/T0ZCU6GXiMdPMjMzOrUVdzN91HGlHdWeZpmeKHxNrA/wAf\ncG8mM7Oe6fPVTRJXAPdFcHqBYZmZNbUiqpu+HxHHSrqmk6drujNd2STGApsBB5UciplZS+qquuni\n/O9ZnTzX9MWPfCvSs4HjI3iz7HjMzFpR1SQREe33Vni1YhkASbsXGlV9HAK8SbqXtpmZ9UItg+nO\nl7RR+4qkA4BvFBfS4pNYFjgN+KIbq83Meq+WO9PtA1wu6ZPAh4BPATsVGtXiOwq4LYJ7yw7EzKyV\n1dS7SdJ6wBTgaeDjEfFG0YF1EUuXLfQSSwFPAeMieKhxkZmZNa8iejd1/IEdRqqeujuf7AOdvKwZ\nHA7c4wRhZrb4qpYk8pxNVUXEjG4PLl0I7AY8HxEb5W3DSPeVXh2YAewXES/l504k/ci/CxwTETd1\ncsyq2VBiCeAJYP8I7uouPjOz/qIp526S9CHgNeDiiiRxJjAnIs6UdDywQkScIGk08CtgC2AV4GZg\n3YiY1+GYXSWJQ4GDI9ixsDdlZtaCipwqvNci4nbgxQ6b9wQm5uWJwN55eS9gUkTMzaWUJ4AxtZ4r\nj4s4Afj24sRsZmYLFJokqhgeEbPz8mxgeF4eCcys2G8mqURRq48BLwG3LHaEZmYG1NYFtjAREZK6\nqu/q9DlJEypW2yKiDfgC8D2PizAzA0ljgbGLe5wyksRsSStHxCxJI4Dn8/ZngVEV+62aty0iIiZU\nrkusAWxIuv+2mVm/ly+e29rXJZ3Sm+OUUd10NWnKDPK/Uyq2j5c0WNKawDrAPTUe8yBgcgRv1TVS\nM7N+rtCShKRJwPbAipKeIU3ncQYwWdIR5C6wABExTdJkYBrwDnBU1ND1SkLAwSxIPGZmVictfz8J\niS2BS4D13B5hZta5puwC2yAHA5c4QZiZ1V9LlyQkBpMat8dE8PdyIzMza179tSSxC/CoE4SZWTFa\nPUl8igV30DMzszpr2eqmPJnfS8AqEbxUdlxmZs2sP1Y3rQXMcoIwMytOKyeJ9YC/lR2EmVlf1upJ\nYnrZQZiZ9WWtnCTWxSUJM7NCtXKScHWTmVnBWj1JuLrJzKxALZkkJIYCy1BlKnEzM6uPlkwSpPaI\n6Z6vycysWK2aJNweYWbWAK2aJNbF7RFmZoVr1SThkoSZWQM4SZiZWVUtOcEfxJvAyhG8UnY8Zmat\noL9N8PeSE4SZWfFaNUm4qsnMrAGcJMzMrKpWTRLu/mpm1gBNlyQk7SzpMUmPSzq+ym4uSZiZNUBT\nJQlJA4EfATsDo4EDJG3Qya5OEnUiaWzZMfQl/jzrx59lc2iqJAGMAZ6IiBkRMRe4DNirk/1mNDSq\nvm1s2QH0MWPLDqAPGVt2ANZ8SWIV4JmK9Zl520IieKdhEZmZ9WPNliRaa2SfmVkfN6jsADp4FhhV\nsT6KVJpYSBp1bfUi6ZSyY+hL/HnWjz/L8jXVtBySBpEapXcEngPuAQ6IiEdLDczMrJ9qqpJERLwj\n6QvAjcBA4BdOEGZm5WmqkoSZmTWXZmu4BkDShZJmS3qoi31+kAfcPShpk0bG12q6+zwljZX0sqT7\n8+NrjY6xVUgaJekWSY9IeljSMVX28/ezBrV8nv5+1k7SEEl3S3pA0jRJ36myX+3fz4hougfwIWAT\n4KEqz+8KXJeXtwTuKjvmZn7U8HmOBa4uO85WeAArAxvn5WVJbWgbdNjH38/6fp7+fvbsM106/zsI\nuAvYrsPzPfp+NmVJIiJuB17sYpc9gYl537uBoZKGNyK2VlTD5wnQ43nm+6OImBURD+Tl14BHgZEd\ndvP3s0Y1fp7g72fNIuKNvDiY1Lb77w679Oj72ZRJogadDbpbtaRY+oIAtslFz+skjS47oFYgaQ1S\nCe3uDk/5+9kLXXye/n72gKQBkh4AZgO3RMS0Drv06PvZVL2beqjjlYVb4HvvL8CoiHhD0i7AFGDd\nkmNqapKWBS4Hjs1XwIvs0mHd388udPN5+vvZAxExD9hY0nuAGyWNjYi2DrvV/P1s1ZJEx0F3q+Zt\n1gsR8Wp7ETUirgeWkDSs5LCalqQlgCuAX0bElE528fezB7r7PP397J2IeBm4Fti8w1M9+n62apK4\nGvgUgKStgJciYna5IbUuScMlKS+PIXWN7liPaUD+nH4BTIuIc6vs5u9njWr5PP39rJ2kFSUNzctL\nATsB93fYrUffz6asbpI0CdgeWFHSM8ApwBIAEfGziLhO0q6SngBeBw4rL9rm193nCewD/I+kd4A3\ngPFlxdoCtgUOAv4qqf0/30nAauDvZy90+3ni72dPjAAmShpAKgRcEhG/l/RZ6N3304PpzMysqlat\nbjIzswZwkjAzs6qcJMzMrConCTMzq8pJwszMqnKSMDOzqpwkrF+RdFKJ5+5s+g6zpuZxEtavSHo1\nIpZrtXNLGhQR79Q7JrPuuCRhfZakqyT9Od/M5sh8A5al8o1rLpG0uqTHJF0k6W+SLpU0TtKdkqZL\n2iIfZ5ikKXkW0j9J2ihv377iRjh/kbRsvkHObZJ+l4/9k/YpJfJrTss3hPmTpJXytvdKulzSPfmx\nTd4+Icd5B2kU7Yqd7WdWqLJvkOGHH0U9gBXyv0sBDwHDgFcrnl8DmAu8nzQr5p9J91WHNOf+VXn5\nh8DX8/IOwP15+Wpg67y8NGnu/rHAm/nYA4CbgE/kfeYBu+Xl7wIn5+VfAdvm5dVI8xgBTADuBZbs\naj8//Cjy0ZRzN5nVybGS9s7LqwLrdLLP3yPiEQBJjwA35+0Pk37oIc0v9HGAiLhF0n9JWg64EzhH\n0qXAlRHxbC403BMRM/IxJwHbkWY5fTsirs3HvI80+RrAR4ANKgocy0lahjR989UR8VYX+y0dC24y\nY1Z3ThLWJ0kaC+wIbBUR/5F0CzCkk13fqlieB7xdsVz5/2OR+fcj4ruSfgfsBtwp6aPtz3V43by8\nPLfDuQZV7LNlRLxN5QtTMqhMAJ3uZ1Ykt0lYX7U88GJOEBsAW+XtcyX19OLoduBAmJ98XoiI1ySt\nHRGPRMSZpGqh9fL+YyStkWfi3B+4o5vj3wQc074i6YM17rdxD9+HWY85SVhfdQMwSNI04HTgT3n7\nz0nTUl9CuuLv2L0vOlmeAGwm6cF8rEPy9mMlPZS3vw1cn7ffC/wImAY8GRFXVTl2+/oxwOa5YfwR\n4LNV4um432e6+QzMFpu7wJrVUS5pHBcRe5Qdi1k9uCRhVl+dlU7MWpZLEmZmVpVLEmZmVpWThJmZ\nVeUkYWZmVTlJmJlZVU4SZmZWlZOEmZlV9f+l5DekVF2HBgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10850af10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# ideal gas parameters for air\n",
    "gamma = 1.4\n",
    "R = 8.314 * ur.J/(ur.degK*ur.mol)\n",
    "MW = 28.966 * ur.grams/ur.mol\n",
    "\n",
    "# ambient conditions\n",
    "Pa = 1.0 * ur.atm\n",
    "Ta = ur.Quantity(20,ur.degC).to(ur.degK)\n",
    "\n",
    "def Qunchoke(P,Pa):\n",
    "    Pr = Pa/P\n",
    "    rho = MW*P/(R*Ta)\n",
    "    b = 2.0*(gamma/(gamma-1))*rho*P\n",
    "    Q = np.sqrt(b*(Pr**(2.0/gamma)-Pr**((1.0+1.0/gamma))))\n",
    "    return Q.to(ur.kg/ur.m**2/ur.sec)\n",
    "\n",
    "P = np.linspace(1.0,3.0,100) * ur.atm\n",
    "Qu = [Qunchoke(p,Pa) for p in P]\n",
    "\n",
    "plt.plot(P,[q.magnitude for q in Qu])\n",
    "plt.title('Unchoked Flow from an Orifice, CD = 1')\n",
    "plt.ylabel(Qu[0].units)\n",
    "plt.xlabel(P.units)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Choked Flow\n",
    "\n",
    "When the flow is choked, flow through the orific depends only on the upstream pressure $P$\n",
    "\n",
    "$$\\dot{m} = C_D A \\sqrt{\\gamma \\rho P \\left(\\frac{2}{\\gamma + 1}\\right)^\\left(\\frac{\\gamma + 1}{\\gamma - 1}\\right)}$$\n",
    "\n",
    "| Parameter | Description | Units |\n",
    "| :-------: | :-: | :-: |\n",
    "| $\\dot{m}$ | mass flow rate | kg/s |\n",
    "| $C_d$ | discharge coefficient | dimensionless|\n",
    "| $A$ | area of discharge | m**2 |\n",
    "| $\\gamma$ | ratio of heat capacities $\\frac{C_p}{C_v}$ | dimensionless |\n",
    "| $\\rho$ | gas density | kg/m**3 |\n",
    "| $P$ | pressure | Pa |\n",
    "| $T$ | temperature | K |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x10877ba50>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEZCAYAAABiu9n+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4VHXa//H3DVhAkCK9t0UFFVmKrg0Suq7uPtss2FbX\nuvZGUSCCDextXfVZG4Kuu+v6032AkABBxIpdUJFAgFBCERDEhty/P87JOsZMMiknM0k+r+uai5nT\nvvdMhnPPOd9m7o6IiEhx6iQ7ABERSV1KEiIiEpeShIiIxKUkISIicSlJiIhIXEoSIiISl5JENWZm\nGWY2rZz7PmFmkyshhrPNbGFV7VfGMo42s8/MbIeZnRRlWanKzGaa2Rkxr28ys01mts7MOoSfjSUz\nRkltShIpzsxOM7PF4X/mdeF/+qPD1RXp5OIV3L9UZtbZzPaEsRc+3o2yzCImAfe5eyN3f7EKy610\nZraPmd1qZqvMbJeZLTOza0rbz92Pd/dp4TE6AlcBB7l7W3dfE342UX8P9jeze8LYd5jZcjO728wO\nCNfnhe/pCzPbamaLzOyCKJKXmV0S/n/62swer+zj10RKEinMzK4C7gZuAloCHYAHgRMLN6loERXc\nP1GNw5NRI3fvU0VlAnQElha3wkJVGEtF/QNIA0YCDYEzgPPN7N7iNo7z/joCW9x9S6SR/jiOvYG5\nwMHAcHdvBPwC2Az0Dzdz4Jfuvn8Y423AaOBvEYS0FpgMPBbBsWsmd9cjBR9AY2AH8NsStpkI/B14\nEvgC+AjoG7P+YCAH2BquOzFm3ePA5PB5I2A+cE/4+iAgC9gCfAL8Pma/A4AXge3AGwT/4RbGia8z\nsAeoU8y6s2P3A44C3gK2AW8CvwiXpwEfxGyXBbwZ83ohcFIxx88Fvgd2hZ/N3uFncROwKFzeNV65\n4TFywve3KPxbvAg0B6aH7/9NoFMJf59/AOvDYy8Aesase4Ig4f8njO91oGuc4wwGvgLaFVk+ANhd\nuF+R9/cl0C1cdm54jF3hZ7KD4CT5o78P0Cz8XqwFPgf+HVPWL4H3wu/SIuDQBL/HfwI2AA1K2GYl\nkF5kWf8w1l4R/f+aDDyezP/j1eWR9AD0iPOHgRHAd8WdYGO2yQhPHiMIrgpuAV4L1+0FLAfGAPXC\nk+0XQI9w/eMEt2MOCE92k8Ll+wFrgLMIrjQPBzYBB4frnw0f9YFeQD7wcpz4Ck9CdYtZdzZhkghP\nTluBUWGZp4QnqaZhOV+F2+wFFITx7Reu2wU0jVP+j04+4QkzjyB51gFaxSs3ZvtlQBdgf2AJ8BmQ\nDtQlSM6PlfD3OTuMcy+CK8J3Y9Y9QfBrul94rKeBZ+Ic5zZgfpx1ecB5cd5fPYLkf064fiCwppi/\nT2GS+D/gGYIfKPWAY8PlfcLPvT/B9+zM8LPdO4Hv8bOUcjIu+neKWb4KuDDOPn8J/3bFPd5LIK6b\nSotLj+Ch202p6wBgs7vvKWW7he4+24Nv/tNA73D5kcB+7n6bu+929/kEv1pPjdm3HcGJ5e/uPiFc\n9ktgpbs/6e573P094Hng92ZWF/gNMMHdv3L3JQQnytJu22wO7zVvDW+hFXUC8Km7Tw/LfJbgCuYk\nd/+K4Jf+QKAvwa/ZRcAx4Xv8zN23llJ+IQeecPePw891WLxyY7Z/3N1XuvsXwCxgmbvPc/fvCa4U\n4t4+c/cn3P1Ld/8OuBHobWaNYo79vLsvDo81nSAhF6c5wa/x4qwP1//k/bn77iLbxv07mVkbgh8b\nF7r79vA7U9iw4HzgYXd/ywNPAd8QfP6laRbGWB7rCH4o/IS7X+zuTeM84n2OPzpEOWOqdeolOwCJ\nawvQ3MzqlJIoCmKe7wL2NbM6QFuCX9yxVoXLIThhnEBw6+HhmG06AUeYWeyJtx7wFMHJqF6R465O\n4L0cUMp7aFvMcWJjXQAMIrhqWUDwa3EgwYkqJ4HyY8XGXlq58OPP92tgY5HXDYsrJPwb3AL8DmhB\n8Isdgs9wRzHH/iresQiu5LrHWdeW4IqkUNG/eaI6AJ+7+/Zi1nUCzjSzS2OW7QW0SeC4W/jx51kW\n7Qiu7KJQneqjkkpXEqnrNYKT4P+UsE1Jv4bWAR2KVF52IrjfXLjvo0AmMNPMGoTLVwMLivwya+Tu\nfyY4Ge0mqFwsFPu8vNaGscWKjXUBwe2y4wiSQmHSGBg+L4vYz6y0ckvatzSjCK5IBrt7Y4JbVlC+\nk1M2QeJuH7vQzI4A2gPzyhljrDVAMzNrXMy61cDNRb4TDd397wnGPjzm+5UQM+tPkCReibP+r0Va\nzcU+PkygCF1JJEhJIkWFv+gmAA+a2a/MrIGZ7WVmI81sSrhZSSecNwiuLK4L9xtEcCvp2dh93f0S\n4FPgJTPbl+C+dA8zOz3cby8z629mB4W3RZ4HMsysvpn1JKi7qOh/uFlhmaeaWT0zO5mg8vw/4fpX\ngQMJ7om/6e5LCa94gJfLWFbsZzazlHKLbl+WE3xDgiT/uZntR3BVES+OErn7XIIWQv8ys55mVtfM\njgSmAX9x99xyxhhbxnqCv8NfzKxJ+Hc/Llz9KHChmQ0IG03tZ2YnmFlD+G+fm3jNSacRJKB/mdmB\nZlbHzA4ws3FmNrJo3GFz2V8S1I1MC29pFhfvhf5Di7mij0Pjvc/ws9uX4Iq4bti0uG7in1TtoySR\nwtz9LoJ27TcQ3OZYDVwM/LtwE356gvZw328JmsqOJLhd8QBwhrsvK2bf8wlu5bwAfEtwr/4Ugl/U\n64FbCVoHAVxCcALcQNBCprSmhPESyH/L96BJ5i+BqwmuVq4haBL5ebh+F/A2sCTmPvurQJ67by56\n4ETjCY8ft9xi4o/7eRfjKYJbV2sJWpa9VoFjAfyWoBJ6NsHtqmnA/7r7pUW2Ky1hl1TmGQSNJT4h\nuBV2GYC7vw2cR/Ad+pyg8v7MmH07EOcXf/g9HBIeM4sfWsU1I2jRVeglM/uC4Ds+FrgT+GMp76U8\nxhP8eBoNnE5wm+/6CMqpMSyo74zo4GZjCf4Qe4APCf7o+xE02+xE0BLjD+6+LWb7cwiavl3m7nMi\nC05EKizsB/EucFh4pSk1TGRJwsw6E9wrPdjdvzGzvxNc3vciaLUz1cxGEzQ3HBPeuphBcEuhHcG9\nzB4JtO4REZGIRHm76QuCS9cGZlYPaEBQmXoSQbNJwn9/HT7/FUE78e/cPY+gjf+ACOMTEZFSRJYk\nwvu6dxLcY1wHbHP3LKCVuxc2/Ssg6NAEQTO5/JhD5BNcUYiISJJEliTMrBtwBUGvzrZAQzM7PXab\nsANYSfe71ExNRCSJouxM1w94NWy5gpk9TzCw1wYza+3uG8JenoWdk9YStJIo1J5i2qubmRKHiEg5\nuHuZm0hHmSQ+AcabWX2CnqlDCMYI+pKgbf2U8N8Xwu1fBGaY2V0Et5l+Fm7/E+V5o1I8M8tw94xk\nx1FT6POsPPosK1d5f2BHliTc/X0zewpYTNAE9h3gEYIRR58zs3MJm8CG2y81s+cIhnbeDVzsUbbP\nFRGRUkU6dpO7TwWmFln8OcFVRXHb38JPe6aKiEiSqMe15CQ7gBomJ9kB1CA5yQ5AIu5xHQUzc9VJ\niIiUTXnPnbqSEBGRuJQkREQkLiUJERGJS0lCRETiUpIQEZG4lCRERCQuJQkREYlLSUJEROJSkhAR\nkbiUJEREJC4lCRERiUtJQkRE4lKSEBGRuJQkREQkLiUJERGJS0lCRETiUpIQEZG4lCRERCQuJQkR\nEYlLSUJEROJSkhARkbiUJEREJC4lCRERiUtJQkSkhrMhY/uWd996lRmIiIikDht8fVdaLJ1Bv5wB\nzC3fMXQlISJSw1j6+APsd6fOot/Dy6n3zf4svrB/uY/l7pUZ248PbnYg8GzMoq7AeOBp4O9AJyAP\n+IO7bwv3GQucA3wPXObuc4oc093dIgtaRKSasrSJ9WmS9ze6Z57MpoPXsa7/uT5n6hwo/7kz0iTx\no4LM6gBrgQHApcBmd59qZqOBpu4+xsx6AjOA/kA7IBvo4e57Yo6jJCEiEsPSMurSaN1ddJtzETtb\nf8Gaoy732XdN/9E25Tx3VmWdxBBgubuvMbOTgIHh8ieBHGAM8CvgGXf/Dsgzs+UESeX1KoxTRKRa\nsLQMo8GmcfSadz179trDpyeO8Zn331WZZVRlkjgFeCZ83srdC8LnBUCr8HlbfpwQ8gmuKEREJIaN\nvOI8erxyO/turc+Kofews/UYn59R6beGqiRJmNnewInA6KLr3N3NrKQ3VjX3w0REqgEbds1JtH/j\nYXquaEHusGls73ihz7/xm6jKi5skzKwvwQnaKOZE7e7vlKGckcDb7r4pfF1gZq3dfYOZtQE2hsvX\nAh1i9msfLisaW0bMyxx3zylDLCIi1Y4NHXMUbd55isPf7UrusJfISzvT503aHnd7s0HAoAqXG6/i\n2sxyCJJDfaAv8EG46jBgsbv/IuFCzJ4FZrn7k+HrqcAWd59iZmOAJkUqrgfwQ8V1d48JUhXXIlKb\n2OBxPWn10XQ6vtKbFYMXsvnA033eTWvKfJyoWjeZ2fPARHf/MHx9CHCju/82wcD2A1YBXdx9R7is\nGfAc0JGfNoEdR9AEdjdwubtnFjmekoSI1HiWPr4NzT99hq7Zx7H66A/YeOgoz75lSbmPF2GSWOru\nPUtbVlWUJESkJrO0CY1ouvIJumf+mg2981j/87M9a8rCCh83wiTxLLCToAOcAacBDd391PIEWlFK\nEiJSE1naxL3YP/9Bus05h+0dN5P/i4s9847nK+34ESaJ+sBFwLHhopeBh9z96zJHWQmUJESkJrG0\nDGO/gpvoOvdqvt3vG1YdN8Zn3ftQpZeT6j2uK4uShIjUFDbyssvovOBm6n1dj5WDp/Jly4wo+jpA\ntFcSxwATgc780GTW3b1rWQurDEoSIlLd2fCrT6bDqw/QaG0Tcoc9wo52V/j8G7+LtMwIk8SnwBXA\nOwSD7gHg7pvLWlhlUJIQkerKhl6XTtvFj9NySXuWD/8n2zqf4/MnfVklZUc4dtM2d59VjphERASw\nIWP70OqDafz89Z6sGJrN20cP8HmTC0rfM/kSuZK4DagLPA/8t+t3GXtcVxpdSYhIdWGDb+hM86Uz\n6JJzJHnHLWZTr9N87s3LkxJLhLebcih+WI60shZWGZQkRCTVWfr4A2i2fBrdskawrt9nrD/8DM++\n7c2kxqTWTSIiyWVpE+vTeNUjdJ9zKlt6bGBt//N8zu0pcbs+yiuJJgStm44LF+UAk9w97sBSUVKS\nEJFUY2kZdWm47g66Zf2ZXc13svqYK3323U8mO65YUY/d9CHB5EAGnAEc5u6/KU+gFaUkISKpIpj0\nZ/NoOs8fDwYrB2X4zAduT3ZcxYkySbzv7r1LW1ZVlCREJBXYiCvPoeMrd9BgS0Nyh97PzjbXRNUR\nrjJE2QT2KzM71t0XhgUdA+wqa0EiIjWBDbv2BNq98QiHLG/F8uFPs73jBVFO+pNsiVxJHA48BTQO\nF20FznL39yOOLV48upIQkSpnQ8YcSZt3n6Lt291ZPmwmW7uf4fMmbU12XImKvHWTmTUGSFaFdUwc\nShIiUmVs8PU9aLlkBp1e/jkr015l88GjfO5Nq5IdV1lFWSdxKzAlZlKgpsDV7n5DuSKtICUJEakK\nlj6+FQcsm07X7HTyf7GEgkNP9+xbk3IHpTJEmSTec/fDiyx71937lLWwyqAkISJRCib9yXuMbpm/\nYeMha1jX74+eNWV+suOqqCgrruuY2b6F80eE80vsXdaCRERSWTDpz9r76TPnXLa338oHp5/mmXf+\nPdlxJVsiSWI6MNfMHiPoJ/FHgopsEZFqL5z050YOmXcdu+t/y8e/vtJn3fdAsuNKFQlVXJvZSGBw\n+DLL3TMjjarkWHS7SUQqhY28/M90XnAr9b7am5XpU/my1cRU7utQEVHebgL4GNjt7llm1sDMGrn7\njrIWJiKSCmz41b+n/WsPcnB+U3KHPsYX7S+JetKf6iqRiuvzgfOAZu7ezcx6EMxxPbjEHSOiKwkR\nKS8bOjqNtosfo9WHHVk+4nm2dj7H50+qFT94Ix2WAxgAvF7YosnMPnT3Q8sVaQUpSYhIWdmQcYfS\n6oMZdHi1F7lD5rPlwNN93uT1yY6rKkV5u+kbd//GzAoLqkcx80uIiKQaS7+hAy0+foZ+845i1XHv\n8tZFPX3uzZ8kO67qJJEkscDMrgcamNlQ4GLgpWjDEhEpP0uf0JSmuU/Rb84JrP/5chZfcJRn3/Z6\nsuOqjhK53VQXOBcYFi7KBP7XkzRbkW43iUg8ljZxHxqvfphumaeztftG8o+4wOfcrh+1VNHMdGbW\nDOiQrMH9whiUJETkRywtw2i4/na6Zl/K1013seqYa3z2PX9LdlypJMqK6wXAiQS3pt4GNgGL3P3K\n8gRaUUoSIhLLjr/kWrrkZGB7jJVpN7Grxa01ta9DRURZcd3Y3b8wsz8BT7n7RDP7sOwhiohUHhtx\n5Vl0ePUuemxsRO7QB9nZ9hqfn/F9suOqaRJJEnXNrA3wB6Bw5NeEs3Q4R/b/Ar3C/f4IfAb8HegE\n5AF/iBlldixwDvA9cJm7z0m0LBGp+WzYtSNp99ajHPJpG3KHP8vyEX/y+Td+ley4aqpEksQkgsrq\nRe7+ppl1IzjJJ+peYKa7/y5sPrsfcD3B8B5TzWw0MAYYY2Y9gZOBnkA7INvMerj7njKUJyI1kA0Z\n24/W7z1Nn7d6kDtsFquO7e3zJm9Jdlw1XZkqrst88GCionfdvWuR5Z8AA929wMxaAznuflB4FbHH\n3aeE280GMtz99Zh9VSchUovY4Ou703LJM3R6uS8rB73O5p6n+dyb8pIdV3UT9dhN5dUF2GRmjwO9\nCSq+rwBauXtBuE0B0Cp83haIbcucT3BFISK1TDjpz9P0yx5M/pFLWXxBn+o86U91FXWSqAf8HLjE\n3d8ys3sIbi39l7u7mZV0OfOTdWaWEfMyx91zKiFWEUkBljZhP5rmPU7fzN+y8ZA1vHPuEM+aOi/Z\ncVU3ZjYIGFTR48RNEmZ2FPBaBTvN5QP57v5W+PqfwFhgg5m1dvcNYaX4xnD9WqBDzP7tw2U/4u4Z\nFYhJRFKQpWXUpVH+fRw+53x2tNvGB6eP8sw7n012XNVV+OM5p/C1mU0sz3Hi1kmY2V+BI4BlwCxg\ntrtvKHMBZi8Df3L3ZeEVQINw1RZ3n2JmY4Am7l5YcT2DYEDBdkA20D02UalOQqRmsbQMo8HG8XSd\nN5bd++xm1cDxPvO+e5IdV00TZWe6g4GRBMNyNAHmAbMJWjuV2ibZzHoTNIHdG8glaAJbF3gO6MhP\nm8COI2gCuxu4vOgER0oSIjWHjbz8QjotvI29v6jPiiF38WXrceoIF42qGpajAZBGkDR+4e59y1pg\nRSlJiFR/Nvya39D+9b/QJK85y4c/wRftL9KkP9GqkiSRCpQkRKovGzr6WNq88wSt3+9M7rCX2Nr1\nLJ83aXuy46oNlCREJGXZkHG9aPnhdDq8ehgrB7/M5h6jfN5NP2mUItFRkhCRlGPpN7Sj+acz6Dr3\nWFYf8z4Fh4zyubcsTXZctZGShIikDEuf0JimuU/RLeuXbDh8Jet/frZn3fZKsuOqzSo9SZhZR2Aq\nQV+FmcDt7v5duO4Fd/91BeItNyUJkdQVTPqz5iG6ZZ7J1q6bgkl/7ngx2XFJNMNyPEbQ+e0Ngpnp\nFpjZSe6+mWD0VhERoHDSnw23cVj2FXzd+CuW/v4in3XPo8mOSyqupCTRwt3/Gj6/xMxOB142sxOr\nIC4RqSbs+Euv4uAFk6jzXR2WD5/Erha3qK9DzVFSkqhnZvu6+9cA7v60mW0gGDZ8vyqJTkRSlo24\nahQdFt3LgQX7kzv0YXa0u0KT/tQ8JdVJXAW8U3TwPDPrA0x196HRh1dsXKqTEEkiG3bdMNq+9b+0\n+Lgdy0c8y7ZOmvSnGlDrJhGJlA0Z25fW70+j3RsHsWJoFlt+drrPm7wp2XFJYpQkRCQSNvj6rrRY\nOoPOCwaQN/BNNvU8zefevCLZcUnZKEmISKWy9PEH0Gz503TLGs7aAZ+wofcZnn3r28mOS8pHSUJE\nKoWlTaxPk7y/0T3zZDYdvI51/c/1OVPnJDsuqZhIkkQ4p8RFwIPufnEF4qs0ShIi0Qgm/Vl7D92y\nLmBn6y9Yc9TlPvuu6cmOSypHpXemM7NOwCvA/wOeM7NO7r6qAjGKSAoKJv3ZdAO95o1jT73v+fTE\nMT7z/ruSHZekhpL6SQwimEr0UOBNgomCnqyCmESkitjIK86jx8Lb2Xd7fVYMuZudrceqI5zEKu12\n06PADcBkdz+/yqIqgW43iVScDb/m17R74yGarmhB7vCn2N7hIp9/4zfJjkuiE1WdRFt3X2dmbdx9\nfYUirCRKEiLlZ0PHHBNM+vNeF036U7uo4lpE4rLB43rS6qPpdHylNysGL2Tzgadp0p/aJYqhwjsB\nxwJ/AJ4DFqZCxbWShEjigkl/lk2na/ZxrD76AzYeOsqzb1mS7Lik6kUxVPggVHEtUi1Z+oTGNFnx\nBP3mnMSGw1bx9nlpnjVlQbLjkupHFdciNYilTdyL/fMfonvm2WzrtJn8I//smXf+K9lxSfKp4lqk\nFrO0DGO/DbfQde6VfNvoG1YdO9pn3fvX0veU2kLDcojUUnb8ZVfQOWcydb+tx4r0KexqeaP6OkhR\nShIitYwNv/oUOiy6n0brmpA77BF2tL9Mk/5IPFHdbqoDHOnur1YkuMqkJCG1nQ29Lp22ix+n5ZL2\n5A7/F1s7/9HnT/oy2XFJaovsSsLM3nP3w8sdWSVTkpDayoaM7UOrD6bR/vWerBiazZafneHzJhck\nOy6pHqJMEncArwP/8hS4N6UkIbWNDb6hM82XzqBLzpHkHbeYTb1O87k3L092XFK9RJkkdgINgO+B\nr8PF7u77JxhYHvBFuP937j7AzJoBfwc6AXnAH9x9W7j9WOCccPvL3H1OkeMpSUitEE76M41uWSNY\n238ZGw4/3bNvXZzsuKR6StmKazNbCfR1989jlk0FNrv7VDMbDTR19zFm1hOYAfQH2gHZQA933xOz\nr5KE1GjBpD+rHqVb5qlsOXAda/uf73Nun5XsuKR6i6LHdeGB6wCjgC7uPsnMOgKt3f3NssRX5PVJ\nwMDw+ZNADjAG+BXwjLt/B+SZ2XJgAMHtLpEaLZj0Z92d9M66mF3Nd/LRKef47Ls1yoEkValJAvgL\nsAdIByYBO8Nl/RIsw4FsM/seeNjdHwVauXthhVsB0Cp83pYfJ4R8gisKkRormPRn82h6zh8PBstO\nuN5nPnB7suMSgcSSxBHu3sfM3gVw98/NbK8ylHG0u683sxZAlpl9ErvS3d3MSrrn9ZN1ZpYR8zLH\n3XPKEI9IyrARV57Dz165gwZbGpI79D52trlWHeGkMpjZIIIx+CokkSTxrZnVjSm4BcGVRUIKh/Nw\n901m9m+C20cFZtba3TeYWRtgY7j5WoJBBQu1D5cVPWZGouWLpCIbdu0JtHvjEQ5Z3orlw59me8cL\nNOmPVKbwx3NO4Wszm1ie49RJYJv7gX8DLc3sFmARcGsiBzezBmbWKHy+HzAM+BB4ETgr3Ows4IXw\n+YvAKWa2t5l1AX5GMAKtSI1gQ8YcaWcMX8bhT7zEF+3f5Z0/tfAXHj9bCUJSVUKtm8zsYGBw+HKu\nu3+c0MGDE/2/w5f1gOnufmvYBPY5oCM/bQI7jqAJ7G7gcnfPLHJMtW6SascGX9+Dlh89Q6eFfViZ\nvojNB53uc29K+vwsUntE2U9imrufUdqyqqIkIdWJpY9vwwGfPk3XuWnk/+IjCg4b5dm3fJjsuKT2\niawJLHBIkYLqAX3LWpBIbWJpExrRNO8x+s3+DQWHreKdPw32rCnzkx2XSFnFTRLhbZ+xQH0z2xGz\n6jvgkagDE6mOgkl/1t5Pnznnsr39Vt4/4xTPvPMfyY5LpLwSud10m7uPqaJ4SqXbTZKKgkl/Cm6k\ny7zr2F3/W/KOG+ez7nsg2XGJFIqyTqIyelxXGiUJSTU28vI/03nBrdT7am9Wpk/ly1YT1ddBUk2U\nSeKvhD2u3f2gsGXSHHdPtMd1pVKSkFRhw6/+Pe1ff4DGa5qRO/Qxvmh/ic+/8btkxyVSnCgrriva\n41qkRrGhowfSdvET9P6wI8tHPM+Kwef4/Ek7St9TpPqJvMe1SE1hQ8YdSqsPZvDz13qxYvB8Fl9w\nlM+bvD7ZcYlEKZEkUbTH9e+AGyKNSiSFWPoNHWjxyQz6zTuaVce9w+ILD/K5Ny9LdlwiVSHSHtdR\nUJ2EVBVLn9CUprlP0W3OCazv+xnrDz/Ls2/TsPVSLUU66ZCZNSUYQqMe4ais7v5OWQurDEoSEjVL\nm7gPjVc/TLfM0/m8ewFrj7jQ59z+UrLjEqmIKFs3TQbOBlYQUxfh7mllLawyKElIVCwtw2i4/na6\nZV3GV82+ZNWxV/vsux9LdlwilSHKJLEMOMTdvy1vcJVJSUKiYMdfci1dcjKwPcbK9Ensaj5FfR2k\nJomyCewSoCnBDHIiNYqNuPIsOi66ix6bGpE79EF2tr3G52d8n+y4RFJFIlcS/YH/B3wEFI557+5+\nUsSxxYtHVxJSYTbs2pG0e+tRDljWmtxhz7Kt03k+/8avkh2XSFSivN30MfAQQZIorJNwd19Q5igr\ngZKEVIQNGTOANu89RZu3e5A7dCZbu5/h8yZtTXZcIlGLMkm85e79yx1ZJVOSkPKwwdd3p+WSZ+j0\ncl9WDnqdzT1P1aQ/UptEmSTuIrjN9CI/3G5SE1ipFix9fCsO+GwaXbOGkH/k0mDSn1vfT3ZcIlUt\nyiSRQ9g3IpaawEoqs7QJ+9E073G6Zf6Wjb3yWdfvXM+amp3suESSJdLOdKlESUJKYmkZdWmUfx/d\n5pzPjrbbWHP0pZ5557PJjksk2ZQkpFaztAyjwcaJdJ03mt377GbVwPE+8757kh2XSKqIsp+ESEqz\nkZdfyIELb2PvHfuyYvDdfNl6nDrCiVSOuFcSZtbW3ddVcTyl0pWEFLLh1/yG9q/9hcarm5M77DG+\naP9nTfqOJbJIAAAUVUlEQVQjUrxKv91kZrOAZsB8YDbwirvvrlCUlUBJQmzo6IG0eftxWn/QieXD\nXmRr1zM16Y9IySKpkzCz+sAgYARwNLAGmAXMdvfV5Qu1YpQkai8bMq4XLT+cQcdFh7JiSA6bDxyl\nSX9EElMlFddm1hUYCQwHWrv7gLIWWFFKErVPOOnPdLrMO4bVx7xPwSGjfO4tS5Mdl0h1UuWtm8xs\nH3f/pvQtK5eSRO1h6RMa0zR3Gt2yTmB9nxVs6HOWZ932arLjEqmO1ARWaoxw0p+/0i3zDLZ220T+\nERf4nDteTHZcItWZkoRUe8GkPxtuo2v2FXzd+CtWH3utz7rn0WTHJVITlPfcWSeKYGKZWV0ze9fM\nXgpfNzOzLDNbZmZzzKxJzLZjzewzM/vEzIZFHZukDjv+0qs4+F876LTgUpYPn8SyE5sqQYgkXyJj\nN50ITAI680PnO3f3/RMqwOwqoC/QyN1PMrOpwGZ3n2pmo4Gm7j7GzHoCM4D+QDsgG+jh7nuKHE9X\nEjWIjbhqFB0W3UvDgv3JHfowO9pdoUl/RCpflAP85QL/A3xU9ISdQFDtgSeAm4Gr3P1EM/sEGOju\nBWbWGshx94PMbCywx92nhPvOBjLc/fUix1SSqAFs2LUjaLv4UVp83JblI55lW6c/adIfkehEOSxH\nPrCkrAkidDdwLRB71dHK3QunQi0AWoXP2wKxCSGf4IpCahAbMrYvrd97mj5vHUju0ExWH/Nznzd5\nU7LjEpHiJZIkRgOzzGw+8G24zN39rpJ2MrNfAhvd/V0zG1TcNu7uZlbSpUyx68wsI+ZljrvnlBSL\nJJ8Nvr4rLZbOoN+CAeQNeoPFF3T3uTevSHZcIjVVeN4dVNHjJJIkJgM7gH2Bvctw7KOAk8zs+HDf\n/c1sGlBgZq3dfYOZtQE2htuvBTrE7N8+XPYT7p5RhjgkiSx9fAuaLZ9Gv6xhrO3/KYsv6O/Zt76d\n7LhEarrwx3NO4Wszm1ie4yRSJ/GRux9SnoPHHGMgcE1YJzEV2OLuU8xsDNCkSMX1AH6ouO7uRQJU\nnUT1YGkT69Mk7290zzyZTT3zWdfvPJ8zdU6y4xKpraKsk5hpZsPdPbMcccUqPNnfBjxnZucCecAf\nANx9qZk9BywFdgMXF00QkvqCSX/W3sPhWRews/UXfHjamT77runJjktEyieRK4mdQAOC+ojCYZgT\nbgJb2XQlkZqCSX823UCXeePYs9f35A0c7zPvvzvZcYlIQD2uJWls5BXn03HhVPbdXp8VQ+5mZ+ux\nmvRHJLVEmiTMrCnwM4IKaADc/eWyFlYZlCRShw2/5te0e+Mhmua2IHfEU2zvcJHPv7HKB30UkdJF\n2ZnuPOAygpZH7wJHAq+5e3p5Aq0oJYnks6FjjqLNO0/R+t2u5A57ia3dzvR5k7YnOy4RiS/KJPER\nwVAZr7n74WZ2EHCru/9P+UKtGCWJ5LHB43rS6qPpdHylNyvTX2HTQaN83k1rkh2XiJQuytZNX7v7\nV2aGme3r7p+Y2YHliFGqKUu/oR3Nl02nf/ZxrD76AxZfeKhn37Ik2XGJSPQSGpYjrJN4Acgys60E\nTVelhrP0CY1psuIJ+s05iQ2983j7vIGeNWVhsuMSkapT1ulLBxGMwzTb3b8tZfNI6HZT9Cxt4l7s\nn/8Q3TPPZlvnzeQfebFn3vF8suMSkfKLpE7CzOoRjP56UEWCq0xKEtGxtAxjvw230HXulXzb6GtW\nHTvGZ93712THJSIVF0mdhLvvNrNPzayTu68qf3iS6uz4y67goAWTqfdNPXKH3sKulpPV10FEEqmT\naAYsMbM3gS/DZe7uJ0UXllQVG371KXRYdD8Hrm9C7rCH2dHuck36IyKFEmkCO6i45ckanlu3myqH\nDb0unXaLH6fFkvYsH/EPtnU61+dP+rL0PUWkOtKwHJIQGzK2D60+mEb713uyYmg2W342SpP+iNR8\nUXam21HM4u3AW8DV7l6lE8coSZSPDb6hM82XzqBLzpHkHbeYTb1O87k3L092XCJSNaJMEjcBa4Bn\nwkWnAN0Ihui40N0HlbXQilCSKBtLH38AzZZPo1vWCNb2X8aG3md69m1vJjsuEalaUSaJD9z9sCLL\n3guH6Hjf3XuXtdCKUJJIjKVNrE/jVY/Qfc6pbOmxnrX9z/c5t89KdlwikhxRDsuxy8xOBv4Rvv4d\n8HX4vHpVaNQCwaQ/6+6kd9bF7Gqxg49OPtdn3/1ksuMSkeopkSuJbsC9BKO/ArwOXEEw/3Rfd38l\n0gh/Go+uJIoRTvozli7zb8DrOCsHTfKZD0xJdlwikhrUuqkWsxFXnkPHV+6gwZaG5A69j51trlVH\nOBGJFWWdRAfgPuCYcNHLwOXunl/mKCuBksQPbNi1J9L+jYdpurwlucOnsb3jhZr0R0SKE2WSyAam\nA0+Hi0YBo9x9aJmjrARKEuGkP63ffZI273Qjd+hMtnY/w+dN2prsuEQkdUWZJH7SgikZrZpiyq61\nScIGX38QLT+aQaeFh2vSHxEpiyhbN20xszOAGYAR9JPYXNaCpPwsfXwbDvj0afrNTSP/Fx+x+MLe\nnn3Lh8mOS0RqvkSuJDoD9/ND66ZXgUvdfXWkkcWPp9ZcSVjahEY0Xfk43TP/h4LDVrGu7x89a8qC\nZMclItWPWjfVIOGkPw/Qbc65bO+4hfwjL/HMO/9R+p4iIsWr9CRhZveXsJ+7+2VlLawy1OQkEUz6\nU3AjXeZdx+7635I3cKzPuvfBZMclItVfFHUSbxP0qC7uoNXr8qMasJGX/5mDFtxKva/2ZmX6HXzZ\narz6OohIsul2U5LZ8Kt/T/vXHqRxflNyh/2NL9pd6vNv/C7ZcYlIzRLF7aZ73f1yM3upmNVJm5mu\npiQJGzo6jbaLH6PlRx3JHf4vtnY+1+dPKm5YdhGRCosiSfR197fjzEzn7l5iKxsz2xdYAOwD7A38\nP3cfa2bNgL8DnYA84A/uvi3cZyxwDvA9cJm7zynmuNU6SdiQcYfS6oMZtH+tFyuGzGNLj1E+b3JB\nsuMSkZotys50fd397SLLfunu/0kgqAbuvsvM6gGvANcAJwGb3X2qmY0Gmrr7GDPrSdAXoz/QDsgG\nerj7niLHrJZJwgbf0InmH0+ny/yjWHXcO2zsdZrPvXlZsuMSkdqhvOfOOgls86iZHRpT0KnAhEQO\n7u67wqd7A3WBrQRJonDo6ieBX4fPfwU84+7fuXsesBwYkEg5qczSJzS13572H/o+vJK9d7Zk8flH\n+TMv9FOCEJHqIJEe178D/mlmpwHHAmcCCY3bZGZ1gHcIZrJ7yN2XmFkrdy+8vVIAtAqftyUYhrxQ\nPsEVRbVkaRP3ofHqR/l55mls6V7Ae2ef6HNu/79kxyUiUhalJgl3XxFePbwArAKGx1whlLbvHuBw\nM2sMZJpZWpH1bmYl3e+qXk2vCCf9abh+Kr2zLmXXATv56OTzffbdjyU7LhGR8oibJMys6NhAzQhu\nT70R3ts6rJjdiuXu283s/4C+QIGZtXb3DWbWBtgYbrYW6BCzW/twWXGxZcS8zHH3nERjiYqlZRj1\nN19Hz5wJ4PDZ8RPY1XyK+jqISDKEjY4GVfg4JbRu6lzSjmG9QfwDmzUHdrv7NjOrD2QCNwLDgS3u\nPsXMxgBNilRcD+CHiuvuXiTAVKy4thFXnk3HV+6iweaG5A59kJ1tr/H5Gd8nOy4RkUKV3uO6tCSQ\ngDbAk2G9RB1gmrvPNbN3gefM7FzCJrBheUvN7DlgKbAbuLhogkg1NuzaE2j35iMcsqw1ucOfYVun\n83z+jV8lOy4RkcqiHtfliWHImAG0eW8abRf/jNxhs/i825k+b/KWZMYkIlISjQJbFWUPvr47LZc8\nQ6eX+7Jy0Ots7nmaz70pLxmxiIiUhZJElGWmj2/FAcuepmv2YPKPXErBYaM8+9b3qzIGEZGKUJKI\noqy0CfvRNO9xumX+lo298lnX71zPmppdFWWLiFQmJYnKLCMtoy6N8u+j25zz2dFuG2uOutQz73w2\nyjJFRKKkJFEZx07LMBpsnEjXeaPZve935B13g8+6774oyhIRqUpRTDpUq9jIyy/iwIW3sfeOfVgx\n5E6+bHWDOsKJSG1X668kbPjVv6X96w/SeHVzcoc9xhft/6xJf0SkptHtprIeZ+jogbR5+3Faf9CJ\n5cNfYGuXszXpj4jUVEoSie4/ZFwvWn44g46LDmXFkBw2HzjK501eX5kxioikGiWJ0vZLv6EDLT6Z\nTpd5x7D6mPcpOPRUn3vzJ1HEKCKSapQk4m2fPqExTXOfolvWL1nfZwUb+pzlWbe9GmWMIiKpRkmi\n6HZpE/eh8ZqH6JZ5Jlu7biL/iAt8zh0vVkWMIiKpRkmicH1ahtFww210zb6Crxt/xepjr/VZ9zxa\nlTGKiKQaJQnAjr/0KjrnTKLO7jqsTL+VXS1uUl8HEZFaniRsxFWj6LDoXhoW7E/usIfY0fYqTfoj\nIvKDWpkkbNh1w2j71t9o/klbcoc/y7ZOf9KkPyIiP1WrkgSDx/Sn9XvTaPfWgeQOzeTz7qdr0h8R\nkfhqV5IY3dTJG/QGmw4e5XNvXpHsmEREUl3tGuBv8QX9PfvWt5MdhohITVctrySSPce1iEh1U95z\nZ50oghERkZpBSUJEROJSkhARkbiUJEREJC4lCRERiUtJQkRE4lKSEBGRuJQkREQkrkiThJl1MLP5\nZrbEzD4ys8vC5c3MLMvMlpnZHDNrErPPWDP7zMw+MbNhUcYnIiIli/pK4jvgSnfvBRwJ/NnMDgbG\nAFnu3gOYG77GzHoCJwM9gRHAX8xMVzsRMrNByY6hJtHnWXn0WaaGSE/A7r7B3d8Ln+8EPgbaAScB\nT4abPQn8Onz+K+AZd//O3fOA5cCAKGMUBiU7gBpmULIDqEEGJTsAqcI6CTPrDPQB3gBauXtBuKoA\naBU+bwvkx+yWT5BUREQkCaokSZhZQ+BfwOXuviN2nQcjDJY0ymD1GoFQRKQGiXyocDPbiyBBTHP3\nF8LFBWbW2t03mFkbYGO4fC3QIWb39uGyosdU4qhEZjYx2THUJPo8K48+y+SLdKhwMzOCOoct7n5l\nzPKp4bIpZjYGaOLuY8KK6xkE9RDtgGygu1e38cxFRGqIqJPEMcDLwAf8cNtoLPAm8BzQEcgD/uDu\n28J9xgHnALsJbk9lRhagiIiUqNpNOiQiIlUnJfsgmNljZlZgZh+WsM19Yae7982sT1XGV92U9nma\n2SAz225m74aPG6o6xuoiXgfRYrbT9zMBiXye+n4mzsz2NbM3zOw9M1tqZrfG2S7x76e7p9wDOJag\nueyHcdYfD8wMnx8BvJ7smFP5kcDnOQh4MdlxVocH0Bo4PHzeEPgUOLjINvp+Vu7nqe9n2T7TBuG/\n9YDXgWOKrC/T9zMlryTcfSGwtYRN/tsZz93fAJqYWasStq/VEvg8ATRveAK8+A6ibYtspu9nghL8\nPEHfz4S5+67w6d5AXeDzIpuU6fuZkkkiAe2ANTGv8wmay0r5OHBUeOk5M2xlJqUo0kE0lr6f5VDC\n56nvZxmYWR0ze4+go/J8d19aZJMyfT8j7ycRoaK/LFQDX37vAB3cfZeZjQReAHokOaaUFnYQ/SdB\nC7ydxW1S5LW+nyUo5fPU97MM3H0PcLiZNQYyzWyQu+cU2Szh72d1vZJIqNOdJMbddxReorr7LGAv\nM2uW5LBSVkwH0af9hw6isfT9LIPSPk99P8vH3bcD/wf0K7KqTN/P6pokXgTOBDCzI4Ft/sNYUFJG\nZtYq7PiImQ0gaBpd9D6m8N8Oon8Dlrr7PXE20/czQYl8nvp+Js7MmhdOvWBm9YGhwLtFNivT9zMl\nbzeZ2TPAQKC5ma0BJgJ7Abj7w+4+08yON7PlwJfAH5MXbeor7fMEfgdcZGa7gV3AKcmKtRo4Gjgd\n+MDMCv/zjSPoGKrvZ9mV+nmi72dZtAGeDKdYqEMwHNJcM7sAyvf9VGc6ERGJq7rebhIRkSqgJCEi\nInEpSYiISFxKEiIiEpeShIiIxKUkISIicSlJSK0STmqVrLKLG75DJKWpn4TUKma2w90bVbeyzaye\nu++u7JhESqMrCamxzOzfZrY4nMzmvHAClvrhxDXTzKyTmX1iZo+b2admNt3MhpnZIjNbZmb9w+M0\nM7MXwlFIXzOzQ8PlA2MmwnnHzBqGE+S8bGb/CY/9UOGQEuE+N4UTwrxmZi3DZS3M7J9m9mb4OCpc\nnhHG+QpBL9rmxW0nEqlkT5Chhx5RPYCm4b/1gQ+BZsCOmPWdge+AXgSjYi4G/hauOwn4d/j8fmB8\n+DwNeDd8/iLwi/B5A4Kx+wcBX4XHrgPMAX4bbrMHOCF8PgW4Pnw+Azg6fN6RYBwjgAzgLWCfkrbT\nQ48oHyk5dpNIJbnczH4dPm8P/KyYbVa6+xIAM1sCZIfLPyI40UMwvtBvANx9vpkdYGaNgEXA3WY2\nHXje3deGFw1vunteeMxngGMIRjn91t3/Lzzm2wSDrwEMAQ6OueBoZGb7EQzf/KK7f1PCdg38h0lm\nRCqdkoTUSGY2CBgMHOnuX5vZfGDfYjb9Jub5HuDbmOex/z9+Mv6+u08xs/8AJwCLzGx44boi++0J\nn39XpKx6Mdsc4e7fErtjkAxiE0Cx24lESXUSUlPtD2wNE8TBwJHh8u/MrKw/jhYCo+C/yWeTu+80\ns27uvsTdpxLcFjow3H6AmXUOR+I8GXillOPPAS4rfGFmvRPc7vAyvg+RMlOSkJpqNlDPzJYCtwCv\nhcsfIRiWehrBL/6izfu8mOcZQF8zez881lnh8svN7MNw+bfArHD5W8ADwFIg193/HefYha8vA/qF\nFeNLgAvixFN0u/NL+QxEKkxNYEUqUXilcbW7n5jsWEQqg64kRCpXcVcnItWWriRERCQuXUmIiEhc\nShIiIhKXkoSIiMSlJCEiInEpSYiISFxKEiIiEtf/B6rwRkIis8fPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x105bc40d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# ideal gas parameters for air\n",
    "gamma = 1.4\n",
    "R = 8.314 * ur.J/(ur.degK*ur.mol)\n",
    "MW = 28.966 * ur.grams/ur.mol\n",
    "\n",
    "# ambient conditions\n",
    "Pa = 1.0 * ur.atm\n",
    "Ta = ur.Quantity(20,ur.degC).to(ur.degK)\n",
    "\n",
    "def Qchoke(P):\n",
    "    rho = MW*P/(R*Ta)\n",
    "    a = gamma*P*rho\n",
    "    Q = np.sqrt(a*(2.0/(gamma+1))**((gamma+1)/(gamma-1)))\n",
    "    return Q.to(ur.kg/ur.m**2/ur.sec)\n",
    "\n",
    "P = np.linspace(1.0,3.0,100) * ur.atm\n",
    "Qc = [Qchoke(p) for p in P]\n",
    "\n",
    "plt.plot(P,[q.magnitude for q in Qc],P,[q.magnitude for q in Qc])\n",
    "plt.title('Choked Flow from an Orifice, CD = 1')\n",
    "plt.ylabel(Qc[0].units)\n",
    "plt.xlabel(P.units)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Critical Pressure Ratio\n",
    "\n",
    "The critical pressure ratio is the ratio of upstream to downstream pressure of a compressible gas at which leads to choked flow.  Given the ratio of heat capacities\n",
    "\n",
    "$$\\gamma = \\frac{C_p}{C_v}$$\n",
    "\n",
    "the critical ratio is\n",
    "\n",
    "$$\\left(\\frac{P_2}{P_1}\\right)_{critical} =  \\left(\\frac{\\gamma+1}{2}\\right)^\\left({\\frac{\\gamma}{\\gamma-1}}\\right)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Species           CR      1/CR  \n",
      "air (dry)         1.8929  0.5283\n",
      "argon             2.0548  0.4867\n",
      "butane            1.7079  0.5855\n",
      "carbon dioxide    1.8324  0.5457\n",
      "chlorine          1.8506  0.5404\n",
      "helium            2.0488  0.4881\n",
      "hydrogen          1.8990  0.5266\n",
      "nitrogen          1.8953  0.5276\n",
      "oxygen            1.8929  0.5283\n",
      "propane           1.7410  0.5744\n"
     ]
    }
   ],
   "source": [
    "gamma = dict()\n",
    "\n",
    "gamma['air (dry)'] = 1.4\n",
    "gamma['argon'] = 1.67\n",
    "gamma['butane'] = 1.096\n",
    "gamma['carbon dioxide'] = 1.30\n",
    "gamma['chlorine'] = 1.33\n",
    "gamma['helium'] = 1.660\n",
    "gamma['hydrogen'] = 1.41\n",
    "gamma['nitrogen'] = 1.404\n",
    "gamma['oxygen'] = 1.400\n",
    "gamma['propane'] = 1.15\n",
    "\n",
    "print \"{0:<15s}   {1:6s}  {2:6s}\".format(\"Species\",\"CR\", \"1/CR\")\n",
    "\n",
    "for s in sorted(gamma.keys()):\n",
    "    r = ((gamma[s]+1)/2)**(gamma[s]/(gamma[s]-1))\n",
    "    print \"{0:<15s}   {1:6.4f}  {2:6.4f}\".format(s, r, 1.0/r)\n",
    "\n",
    "\n",
    "def isChoked(Phi, Plo = 101325.0, gamma = 1.4):\n",
    "    r = ((gamma+1)/2)**(gamma/(gamma-1))\n",
    "    return Phi > r*Plo\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Procedure to Calculate the Airbag Vent Mass Flowrate\n",
    "\n",
    "Given the initial airbag pressure $P_i$ and absolute temperature $T_i$, the ambient pressure $P_a$, and gas constant $\\gamma$, and orifice vent area $A_o$:\n",
    "\n",
    "**Step 1.** Calculate the intial air density as\n",
    "\n",
    "$$\\rho_i = \\frac{M P_i}{R T_i}$$\n",
    "\n",
    "For adiabatic (isentropic) compression, subsequent values of the airbag pressure, temperature and air density satisfy\n",
    "\n",
    "$$ \\frac{P}{P_i} = \\left(\\frac{\\rho}{\\rho_i}\\right)^\\gamma =\n",
    "\\left(\\frac{T}{T_i}\\right)^\\frac{\\gamma}{\\gamma-1}$$\n",
    "\n",
    "**Step 2.** Given the ambient pressure $P_a$, compute the discharge coefficient\n",
    "\n",
    "$$C_D = 0.9 - 0.3\\frac{P_a}{P}$$\n",
    "\n",
    "**Step 3.** Compute the critical pressure ratio\n",
    "\n",
    "$$r = \\left(\\frac{\\gamma+1}{2}\\right)^\\left({\\frac{\\gamma}{\\gamma-1}}\\right)$$\n",
    "\n",
    "**Step 4.** Compute the mass flow for either unchoked or choked flow conditions\n",
    "\n",
    "$$ \\dot{m} = \\begin{cases}\n",
    "C_D A_o \\sqrt{\\frac{2\\gamma\\rho P}{\\gamma -1}\\left[\\left(\\frac{P_a}{P}\\right)^\\frac{2}{\\gamma} - \\left(\\frac{P_a}{P}\\right)^{1+\\frac{1}{\\gamma}}\\right]}\n",
    "& \\mbox{if } P \\leq r P_a \\\\\n",
    "C_D A_o \\sqrt{\\gamma \\rho P \\left(\\frac{2}{\\gamma + 1}\\right)^\\left(\\frac{\\gamma + 1}{\\gamma - 1}\\right)}\n",
    "& \\mbox{otherwise}\n",
    "\\end{cases}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x10c610f50>"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEZCAYAAABvpam5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcHFW9/vHPA2HfA8q+KAQVRURkEVDiBrgB4k8WFQFx\nRQVFVMB7JV5Fxe2qCOh1BRQQFxAxsolRQCGigKwiShSCBEFkXxLy/P44Z0gzmcn0TKanumee9+vV\nr1RXV3d9uydV3zpLnSPbREREDMcSTQcQERG9J8kjIiKGLckjIiKGLckjIiKGLckjIiKGLckjIiKG\nLckjRpWkGZIO6vA+PinpX5Ju7+R+ondIulbSi+uyJH1H0r8lXSZpR0k3Nh3jeJPk0UMkzZL0qKTV\n+62/UtJ8SRuMURzTJM2VdH/L4/D6suujU/veADgMeKbtdTq1n7Eg6TmSzquJcP4Ar0+WdKakB+rf\nft9+r79M0o2SHpR00Vj9/eu+vzPAug0lHd1v3ecl3STpPkk3SNpviM9dVdKJkv5Zv9efJB0wVDy2\nn2P7N/XpjsDLgXVsb2f7EtvPbP/bRTuSPHqLgb8BT5xEJG0OLEcHT9iDxHGa7ZVaHp8fo31vANxt\n++6BXpQ0aYziGA2PAacDg5XUjgceAZ4KvAk4UdJmAJLWAH4MfBRYDbgC+EEng5W0vqQvSlq+Pn+O\npM9K2lbSUcCkuv5Fko6sb3sAeI3tlYH9gS9LeuEgn780cCGwPrAdsDLwIeAzkj4wyHsG+ntvCMyy\n/ciIv2wMzXYePfIAbqGcLGa2rPs8cBQwH9igrns1cCVwL/AP4OiW7ZcFvgfcBdwDzASeWl87APgr\ncB8lSb1xkDimAacM8tqvgLfWZQH/BcwC5gAnASvX104CDqvL69b4D67PN6YkiP6f/XLgIeBx4H7g\n25QTxXzgrcDfgRlD7Hejuv0B9be5G3gXsDXwp/qbHLeIv8E2wO/qdrcDxwFLtbw+H3gncFPd5qtt\n/F03Aeb3W7cC8CiwScu6k4BP1+V3AJe0vLZ8/W02HWQfk4HvALOBfwNn1vVTgdsoJ+k763faA3hV\n/Q53A0e0fM4OlKQ1q373yXX9bsAFwA3AJ4HlBonjp31/9wFeO6j+vZbrt36v+vdesT6fBXy4/r0e\nBpas615WP+NhYF59z9H1O97a8nnrAz+p3/eu1r93/X90ff2NzqUeU3ks/EjJo/dcBqws6ZmSlgT2\npiSDVg8Ab7a9CiWRvFvS7vW1/SlXdOtRTijvBB6WtALwZWBXl6vEFwJXLWasB9b9TQWeDqwIfLW+\nNqOuB9iJkqxe3PK8rwriCbYvBF4J3O5S2nkrJVFQ3/tMYNch9ttnG8pJex/K9z4KeCnwbGCvvvrz\nAcwDDgVWp/xGLwMO7rfNq4EXAM+tn7XLIJ+1KJsC82zf3LLu6hof9d+r+16w/RBwM/CcQT7vFMqF\nw2aUkswXW15bE1gGWBv4GPBNSklnS+BFwMckbVi3FQtKufNZ8PuziPXljdJylCR97SAxvgKYbvvh\nfut/UmNvLbHsQ/m/sKrtx+u+bftblIuB39X/Ix/vF8OSwDmUC7ENKRcup9fXdgeOBF4HrAFcDJw2\nSKwTXpJHbzoFeAvlYLuecjX5BNu/tn1dXb6GcnDsVF9+jHLim+LiStv319fmA5tLWs72HNvXLyKG\nvSTdUx//lrTWANu8CfiC7Vm2H6QcmPtIWoKSHHaUJMoJ6rOUq1pqrL8eZL8LnZSqabYfdqmqWNR+\n+3zC9mO2L6BcoZ5q+y7bt1NOGlsOtBPbf7Q90/Z8238H/o8Fv22fz9i+z/atlJLY8waJeVFWpJQA\nW90PrLSI1++r659E0tqUpPou2/fanmf74pZN5gLH1JPwDygXFV+y/WD9P3A98DxJ6wN7AvvV7/VN\n4COStgU2B95N+b92ISXB9vc14Crb5w/ynVcH/tl/pe15lBLCGn2rgK/Ynm370QE+Z7D/I1AuGtYG\nPlT/vzxq+9L62rsoJbs/254PfLrle0c/SR69x5Tk8SbK1fXJ9DtYah30ryTdKek/lNJFXyP7KcB5\nwOmSZks6VtKkepLdm3IA3S7pHEnPWEQcP7C9Wn1Mtn3HANusTalK6vMPSr34mrb/CjxIObG+iHI1\neLukTSmliMGSx2BubWe/LevmtCw/PMDzhU7CAJI2rb/NPyXdCxzDgt+2T+tv8dBgnzWEByglxFar\nsCBh3D/I6/ezsPWBf9u+d5B93W27r9TQd9Xf//dYwfattg/rKxnYvsb2h21fbvsYSqkM27+x/enW\nHUj6HKXUs9cgMUBJEAt1gqjtGmvU1/vc2n+7Nq0P/L0mh/42pLTJ3CPpHkqVHZTSSfST5NGDbP+D\nUs3zSkqRvr9TgbOA9WyvSrniW6K+d57t/7H9bGB74DWUUgy2z7e9M7AWcCPwjcFCYNFXd31up7Qx\n9NmAcoLpOzH9GngDpc3g9vr8AEoD8HCrzFo7DAy13+F+XqsTKVfim9RqwY/SmePoJmCSpE1a1m0B\nXFeXr6vPAajVjhu3vN7qVmCypFVGKzjbBw6w7u/9q4lqbB8HdgF2tv3AIj72QuCVfQ3yLV5Paf+5\nrHV3w48aKL/FBrX6qr9/AO9ouShazfYKti8bYNsJL8mjdx0EvHSA+mEoV7r32H5M0jbAG6kHm6Sp\nkjavB8/9lCqLxyU9VdLu9SQ0l1IqeHyQfbeTOKDUF39A0kaSVgQ+BZzectX3a+C9LGjfmFGfX9xy\nJTwSQ+23HYN9xxUpv9tDkp5JqaoZyeeUF6VlgaXr8jKSlgGoJcGfAP8jaXlJOwKvpZQcAc4EniNp\nz/oZR1OqhG7qvw/b/wR+AZxQu8IutYg2nVFVe13tC7zC9j1DbH4KpfH+h7Xb71K1vejLlE4fA5Wq\nhmsmpWrsM/V3XVbS9vW1rwFHtfRoW0XSG0Zhn+NSkkePsv03239sXdWyfDDlpHMf8N88uQvnWsAP\nKT2xrqecsE+h/F/4AKX95G5KVdJgJ8Z27+X4dv3s31BKSg8B72t5/TeUk3Ff8riU0u14ocbyAfa/\nqOdD7bed2Afb5nBKMr6P0t5xer9tB4ptwM+StFGN7dq6zcOU3kp9Dqb8HndSOkW8y/YNALbvolyR\nH0PpGfQCSiPyYPajXBTcSCmBHTJEzKPlGEpV0c1acE/QEQNtaPsxSo+6W4HLKf9HPw8cZfsLw9jn\nQL+56z4epyThTSgljVupVWm2zwKOpVTp3gtcQykxxQC0eBd4bexAmkU50B4H5treRtJkygltQ0oX\nu71s/6dufySlu9zjwCF9jWuStgK+S+l1Md32QA1yERExBsai5GFgqu0tbW9T1x0BXGB7U+CX9Tm1\nuLg3pWFtV0oxu6/YfyJwkO0pwBRJu45B7BERMYCxqrbqX++7G+WGJ+q/e9Tl3Sl3Ls+1PYvSb33b\n2tVwJdsz63Ynt7wnIiLG2FiVPC6UdIWkt9d1a9ru6/kyhwVdKNehNJj1uY3STa7/+tmk+1xERGPG\nYhygHWz/U9JTgAvUb3RL25Y0luMyRUTEYup48qjdBLH9L0lnUu7wnCNpLdt31CqpO+vmsyk9M/qs\nRylxzK7LreufdFc1QJJQRMTw2W63+/2T3tSxB2WwtpXq8gqUrpg7U4ai+EhdfwRlOAcoDeVXUfq9\nP40ySF9fj7DLgW0p7SfTKWMw9d+fO/l9Ruk3mdZ0DIkzcSbOxNgSp0fyvk6XPNYEzqwdpiYB37d9\nvqQrgDNUJg2axYJ+1tdLOoNy/8E8yiirfaWJgylddZejdNU9t8OxR0TEIDqaPGzfwgCDwtn+N+Vm\noIHe8ynKHcH91/+BMvhaREQ0LHeYj70ZTQfQphlNB9CmGU0H0KYZTQfQphlNB9CmGU0H0IYZTQfQ\nSR2/w3wsSbJH0vATETFBjfS8mZJHREQMW5JHREQMW5JHREQMW5JHREQMW5JHREQMW5JHREQMW5JH\nREQMW5JHRESPkLrnnN01gURExMAk1pf4AnCFtNDkeo1I8oiI6FISm0ucDFxdV+1u0xXDgiR5RER0\nEQlJvERiOnAecAOwsc0HbW5tOLwnjMVMghERMQSJScCewIeAlYDPAa+zebTRwAaR5BER0SCJ5YED\ngQ8CtwOfAM6xmd9oYENI8oiIaIDEU4D3UCa6+y2wn82lzUbVvrR5RESMIYmNJY4HbgLWBV5ss0cv\nJQ5I8oiIGBMSL5D4AXA5cC+wmc3bbW5sOLQRSbVVRESH1HsydgE+DGwC/C/wNpv7Gw1sFCR5RESM\nMomlgX2AwwFTek79wGZuo4GNoiSPiIhRIrEy8Hbg/cCfKd1uz++WG/tGU5JHRMRiklgbOISSOC4E\n9rD5Q7NRdVYazCMiRkjimRLfBK6n3Ni3jc0+4z1xQEoeERHDUhvBd6BUSW0HHA9Msbmr0cDGWJJH\nREQbJJYEdqMkjacCXwD2tXmo0cAakuQREbEIEssCb6EMH3Iv8FngTJvHGw2sYUkeEREDkFgNeDfw\nPuCPwDuA34zHnlMjkeQREdFCYgPgA8D+wNnAK2yubTaq7pPeVhERgMQWEt8DrgQeB55rc0ASx8CS\nPCJiwqoTL71M4lxgOvAn4Ok2h9vc1nB4XS3VVhEx4dSJl/4fZcypZSnDh+zerRMvdaMkj4iYMCRW\nAN4KHAbcChwN/LzbJ17qRkkeETHuSTwVeC/wLuASyv0ZlzUbVW9Lm0dEjFsSm0icANxIubFvR5s9\nkzgWX5JHRIw7EttI/BD4HXA38Cybd9nc1HBo40aqrSJiXJBYAnglZfiQjYAvAgfaPNBkXONVx0se\nkpaUdKWkn9XnkyVdIOkmSedLWrVl2yMl/UXSjZJ2blm/laRr6mtf7nTMEdE7JJaWOIDSzfaTwNeB\nTWy+ksTROWNRbXUoZbjivlv6jwAusL0p8Mv6HEmbAXsDmwG7AidIUn3PicBBtqcAUyTtOgZxR0QX\nk1hZ4nDgb8AbKXeFP9/mNJt5zUY3/nU0eUhaD3gV8E2gLxHsBpxUl08C9qjLuwOn2Z5rexZwM7Ct\npLWBlWzPrNud3PKeiJhgJNaROBa4BdgKeK3NzjYXZNypsdPpksf/UuofW/tQr2l7Tl2eA6xZl9eB\nJ93ReRuw7gDrZ9f1ETGBSGwm8W3gWsqNfVvZ7GtzZcOhTUgdazCX9BrgTttXSpo60Da2LWlUrxQk\nTWt5OsP2jNH8/IgYO3XipRdR7gR/AfBVSnvGvxsNrIfV8/HUxf2cTva22h7YTdKrKFcJK0s6BZgj\naS3bd9QqqTvr9rOB9Vvevx6lxDG7Lreunz3YTm1PG72vEBFNqBMv7UFJGqtRJl56g83DjQY2DtQL\n6hl9zyUdPZLP6Vi1le2jbK9v+2nAPsBFtvejDHG8f91sf+Csunw2sI+kpSU9DZgCzLR9B3CfpG1r\nA/p+Le+JiHFEYjmJd1Fu6vsQcCzlHo2vJ3F0l7G8z6OveuozwBmSDgJmAXsB2L5e0hmUnlnzgINt\n973nYOC7wHLAdNvnjmHcEdFhEpMpx/l7gd8DBwEXpwG8e2nB+bn3SbJtDb1lRHQDiY0oXWz3A34K\nfM7m+kaDmmBGet7M8CQRMeYktpQ4FbgCeATY3ObAJI7ekeQREWOiTry0s8QFwM8o84I/3eYj9uCd\nYKI7ZWyriOgoiaWAN1B6Tk0CPg+cavNYo4HFYknyiIiOkFiR0vB9GKVzzEeBX2TipfEhySMiRpXE\nmsD7gHdS7id4g83MRb4pek7aPCJiVEhMkfgacAMwGXihncQxXqXkERGLRWI7SnvGjpQRsJ9h869m\no4pOS/KIiGGrEy+9mpI01qMMH7KfzYONBhZjJskjItomsQzwJsrQIQ8BnwN+lPkzJp4kj4gYksSq\nlAbwQyhDor8XuCjDh0xcSR4RMSiJ9YD3AwcC04FX2VzdbFTRDdLbKiIWIrG5xEmUecGXALa02S+J\nI/okeUQE8MTwIVMlpgPnA38GNrY5zOYfDYcXXSbVVhETXJ14aU9KI/gqlEbwPW0eaTSw6GpJHhET\nlMRywAHABykzen4a+GmGD4l2JHlETDASqwPvqY/LgP1tLm02qug1afOImCAkniZxHPAXYANgqs3u\nSRwxEkkeEeOcxFYSp1Omd70feLbN22xuaDi06GGptooYhyQE7EIZPmQK8CXg7Tb3NxpYjBtJHhHj\nSJ14aR/gcErNwmeB023mNhpYjDtJHhHjgMRKwNuADwA3Ax8BzsvwIdEpSR4RPUxiLcp4U+8Afgm8\n3ub3zUYVE0EazCN6kMQzJb4BXA+sDGxjs3cSR4yVlDwieojE9pRG8O2B44FNbe5qNqqYiJI8Irpc\nnXjpNZSksQ5l4qU32jzUaGAxoSV5RHQpiaVZMPHSw8CxwE8y8VJ0gySPiC5Te069ndJz6kZKg/gv\n03MqukmSR0SXkHgK8D7g3cBFwB42f2g2qoiBDZo8JB3X8tSAWp/bPqRjUUVMIBIbUEa23Q/4EbC9\nzV+ajSpi0RbVVfcP9bEM8HzgJsqAas8Dlu58aBHjW+1u+x3gKuAx4Dk270jiiF4ge9HVqJIuB3a0\nPbc+Xwq4xPa2YxDfsEiybQ29ZURzJLYEjgJ2Ao4Djrf5d7NRxUQ10vNmOzcJrkq5CanPSnVdRAyD\nxA51itdzgN8CT7f5RBJH9KJ2Gsw/A/xR0oz6fCdgWqcCihhP6ui2LwP+izKHxrHA62webTSwiMU0\nZLUVgKS1gW0pDeeX276j04GNRKqtolvUpPEq4L8p84J/Cjgt92hEtxnpebPd5LEusBGlpGIA278Z\n7s46LckjmlbvBt+DUtJYEvgk5ca+xxsNLGIQIz1vDlltJelYYG/KAGytB8Aik4ekZYFfU3prLQ38\n1PaRkiYDPwA2BGYBe9n+T33PkcBb634OsX1+Xb8V8F1gWWC67UPb/4oRnSexJPD/KEnjEUrV7jk2\n85uMK6JT2ultdROwue1h19FKWt72Q5ImAZdQJqjZDbjL9mclfQRYzfYRkjYDTgW2BtYFLgSm2Lak\nmcB7bc+UNB34iu1zB9hfSh4xpmrS2JuSNO4F/gc4N3eDR6/oZG+rvzLC+zps9w3ctjSlCH8PJXmc\nVNefRCniA+wOnGZ7ru1ZlAlttq3tLSvZnlm3O7nlPRGNkJgk8WZKifw9wPspN/f9IokjJoJ2els9\nDFwl6ZfwRA+Rtu4wl7QE8EdgY+BE29dJWtP2nLrJHGDNurwOcFnL22+jlEDm1uU+s+v6iDEnMYky\nzet/A3cCBwMXJWHERNNO8ji7PvoODrUsL5Lt+cDzJK0CnCfpJf1et6QcdNH1avXUviRpRABtJA/b\n35W0DLBpXXVj393m7bJ9r6SfA1sBcyStZfuOWiV1Z91sNrB+y9vWo5Q4Ztfl1vWzB9uXpGktT2fY\nnjGcWCNa1aSxF3A08C/KoIW/StKIXiVpKjB1sT+njQbzqZS2ib/XVRsA+9v+9RDvWwOYZ/s/kpYD\nzgM+DuwC3G37WElHAKv2azDfhgUN5pvU0snllGGpZwI/Jw3m0WG1y+3rKb2m7qUkjwuTNGK86VhX\nXeCLwM62/1x3tClwOmWwxEVZGziptnssAZxi+5eSrgTOkHQQtasugO3rJZ1BaYCcBxzsBZntYEpX\n3eUoXXUXShwRo6He3LcbpdfUo5TRbs9L0oh4snZKHn+y/dyh1nWDlDxipGrS2AX4BOWi6mOU+zSS\nNGJc62TJ4w+Svgl8j9JY/ibgiuHuKKJbSexEuRN8dUrS+Elu7otYtHZKHstS+rHvUFddDJwwkpsG\nOy0ljxgOia2BYyhdyacBp2YYkZhoOja2laQVgEdsP16fLwks03IDYNdI8oh2SDybUtLYmtK28R2b\nYfUgjBgvOnmH+UWUhuo+y1N6QkX0FImnSZxM+T99MTDF5v+SOCKGr53ksYztB/qe2L6fkkAieoLE\nmhJfpbTV/Y2SNL5o83DDoUX0rHaSx4N1VFsAJL0ActBF95NYReITlO7fjwHPtJlmc1/DoUX0vHZ6\nW70f+KGk2+vztSmjiEZ0JYllKPcGHQFMB55vP3GTa0SMgnYng1oKeAalq+6whycZK2kwn9jqUCJv\npNyrcTVwlM11zUYV0d06ORnUCsBhwAa23y5piqRn2D5nJIFGjLaWG/yOBR4E3mxzSbNRRYxv7bR5\nfIdSX7x9fX47pW98ROMkng9cAHyJcq/GDkkcEZ3XTvLY2PaxlASC7Qc7G1LE0CQ2lDiFMlDmj4Dn\n2JyZ4UQixkY7yePROiouAJI2ZsGkUBFjSmJViWMpk4z9DdjU5ms28xoOLWJCaae31TTgXGA9SadS\nhik5oIMxRSxEYingXZS5wn8GbG5z+6LfFRGd0m5vqzWA7erTy2zf1dGoRii9rcaf2hj+WuBzlDll\nDrf5U7NRRYwfHRueRNKOlLGtzgFWA46StOEIYowYFonnAb8EPg0carNzEkdEd2inzeNE4CFJW1C6\n7P4VOLmjUcWEJrGWxLco1aU/BLawyQRgEV2kneQxz/Z8YA/geNvHAyt1NqyYiCSWlTgSuBa4G3iG\nzYlpDI/oPu00mN8v6SjgzcCL6pDsS3U2rJhIarvG6yntGlcB29nc3GxUEbEo7SSPvYF9gbfavkPS\nBpSDPGKxSWwBfJnSnnaQzUUNhxQRbWirt1WvSG+r3iGxBmUMqj0p3cG/keqpiLHXycmgIkaNxCSJ\n91KGSZ8LPCvtGhG9p51qq4hRIbETcBxwF/BSm2sbDikiRmjQkoek/5P0OknpWRWLRWJdiVMpXbw/\nAbwsiSOity2q2urbwPOA6ZIukvSReq9HRFsklpb4EGVujb9Rqqh+mMELI3rfcIYn2RnYFXgucCXw\nC9tndDa84UmDefeQeClwPHAL5e7wvzQcUkQMYKTnzWH3tpIkYCtgF9tdNa9HkkfzJNYBvgC8EDgU\nODsljYjuNWbJo5sleTRHYhLwXsqot18HjrF5qNmoImIoHZuGNmIoEtsBX6MMKbKDzZ8bDikiOizJ\nI0ZMYjXKiLe7AYcDp6WKKmJiGDJ5SJoEvBrYqGV72/5iB+OKLlbHotoX+DxwFrCZzX+ajSoixlI7\nJY+fAQ8D1wDzOxtOdDuJjYETgLWA19lc3nBIEdGAdpLHuraf2/FIoqvVaWAPAz4EHAt8yWZus1FF\nRFPaSR7nS9rF9nkdjya6ksQ2wDeAfwJb29zScEgR0bB2ksdvgTMlLQFPXGna9sqdCyu6gcSKlOFE\n9qWUOtIgHhFAe6PqfhHYDlje9kr1kcQxzknsTGnnWh14js2pSRwR0aedksc/gOvqVLQxzklMplww\nTAXeaZPqyohYSDslj1uAX0k6UtIH6+Owdj5c0vqSfiXpOknXSjqkrp8s6QJJN0k6X9KqLe85UtJf\nJN0oaeeW9VtJuqa+9uXhftEYmsSelPnD7wc2T+KIiMG0mzwuApYGVgRWqo92zAU+YPvZlKqv90h6\nFnAEcIHtTYFf1udI2owy7e1mlEEYT6hjaQGcCBxkewowRdKubcYQQ5B4qsQZlBv+9rJ5n839TccV\nEd1ryGor29NG+uG27wDuqMsPSLoBWJdyR/JOdbOTgBmUBLI7cJrtucAsSTcD20r6O7CS7Zn1PScD\newDnjjS2KCT2Ar5C+U33t3m44ZAioge0c4f5U4EPU0oDy9XVtv3S4exI0kbAlsDlwJq259SX5gBr\n1uV1gMta3nYbJdnMrct9Ztf1MUIST6UMmf4cYA/7Sb97RMQitVNt9X3gRuDpwDRgFnDFcHYiaUXg\nx8Chtp9UHeIyrG968YwhidezYIKmLZM4ImK42ulttbrtb0o6xPavgV9Lajt5SFqKkjhOsX1WXT1H\n0lq275C0NnBnXT8bWL/l7etRShyz63Lr+tmD7G9ay9MZtme0G+t4V3tSHQe8gDK0SJJGxAQjaSql\nN+ViaSd5PFb/vUPSa4DbgdXa+fDa2P0t4HrbX2p56Wxgf8owF/tTBtfrW3+qpC9SqqWmADNtW9J9\nkrYFZgL7UerpF7I4bTTjmcQrKXeJ/4hS2shcGxETUL2gntH3XNLRI/mcISeDqgnjEkqJ4DhgZWCa\n7bOH/HBpR+A3wJ9YUDV1JCUBnAFsQKkG28v2f+p7jgLeCsyjVHOdV9dvBXyX0u4y3fYhA+wvk0H1\nU+8S/xzwSuBAm181HFJEdJGOzCQoaUnKCbwnhl9P8niyOknTKcCllHnE7204pIjoMiM9by6ywdz2\n45RxjaKHSCwl8XHgp8BHbA5I4oiI0dROm8clkr4K/AB4EBClk9QfOxpZjIjEFOB7wD3A82z+2XBI\nETEOtdPmMYMButLafkmHYhqxiVxtVWf3O5DSCeF/gK9mIMOIGEpH2jx6zURNHnUu8a8DzwL2tbm2\n4ZAiokeM9LzZzh3mH2Thkse9wB9sXzXcHcboktiRciPnT4G32DzScEgRMQG0U211KuWmsp9R2jte\nTZnnYUPgR7aP7XSQ7ZpIJQ+JJYGPAgcDb7M5p+GQIqIHdazaStLFwCttP1CfrwhMp4x6+wfbzxpB\nvB0xUZKHxDqU0oaB/eyB77aPiBhKR7rqVk9hwV3mUAYpXNP2Q5AqkrEmsQvwB8ow+a9I4oiIJrTT\nVff7wOWSzqJUW72WMoTICsD1nQwuFpCYROlF9RZgH5tfNxxSRExgbfW2krQ1sAOlmuRS28MaVXes\njNdqq1pNdTrwMPBmm381HFJEjBPpqsv4TB4SL6dM1HQC8CmbzCUfEaOmY111oxkSSwBHUXpTvdnm\nooZDioh4QpJHF6o3/Z0MTAa2TqN4RHSbdnpbxRiS2IIyU+PfgJckcUREN0rJo4tI7EuZ5OpQm1Ob\njiciYjBJHl2gdsP9DLAn8HKbqxsOKSJikZI8GlbnFf8BpRv01jZ3NxxSRMSQ0ubRIIlnU6bkvRp4\nVRJHRPSKlDwaIvFa4FvA4TYnNx1PRMRwJHmMsTpp0+HA+4HX2lzecEgREcOW5DGGJJYGTgS2Al5o\n84+GQ4qIGJEkjzFSG8Z/QplIa0ebBxoOKSJixNJgPgYkng78ljKU+p5JHBHR65I8OkxiO+BS4Dib\nD9o83nRmp734AAAL90lEQVRMERGLK9VWHSSxO/BN4MBMExsR40mSR4dIvBM4mnL/xu+bjiciYjQl\neYyy2hX348AbgRfb3NxwSBERoy7JYxRJLAkcD7wA2N7mzoZDiojoiCSPUVLv4TgFeAplKPX7Gw4p\nIqJj0ttqFEisAJwNLE1p40jiiIhxLcljMUmsDPwCuAN4g80jDYcUEdFxSR6LoU4XewFwLfBWm3kN\nhxQRMSaSPEZI4inARZQbAN9jM7/hkCIixkySxwjUxPFLSnXVB23ccEgREWMqyWOYJNYALgR+Bnw0\niSMiJqIkj2GQWJ2SOKYD/5XEERETVUeTh6RvS5oj6ZqWdZMlXSDpJknnS1q15bUjJf1F0o2Sdm5Z\nv5Wka+prX+5kzIORWAU4HzgPOCqJIyImsk6XPL4D7Npv3RHABbY3pbQbHAEgaTNgb2Cz+p4TJKm+\n50TgINtTgCmS+n9mR9X7OH5OaRw/IokjIia6jiYP2xcD9/RbvRtwUl0+CdijLu8OnGZ7ru1ZwM3A\ntpLWBlayPbNud3LLezpOYhnKJE43A+9P4oiIaKbNY03bc+ryHGDNurwOcFvLdrcB6w6wfnZd33F1\nrKrvAw8Ab0t33IiIotGxrWxb0qheyUua1vJ0hu0ZI/scBHwZmAy8MjcARsR4IGkqMHVxP6eJ5DFH\n0lq276hVUn0jz84G1m/Zbj1KiWN2XW5dP3uwD7c9bZTi/DDwIsqw6o+O0mdGRDSqXlDP6Hsu6eiR\nfE4T1VZnA/vX5f2Bs1rW7yNpaUlPA6YAM23fAdwnadvagL5fy3s6QmI/4GDKIIf3dnJfERG9qKMl\nD0mnATsBa0i6FfgY8BngDEkHAbOAvQBsXy/pDOB6YB5wsO2+Kq2Dge8CywHTbZ/buZh5MfB5yrDq\ng5ZwIiImMi04P/c+Sbatobcc7P08Dfgt8BabC0YvsoiI7jTS82buMK/q0Oo/A45J4oiIWLSUPHii\nS+5ZlIb4d+dejoiYKEZ63sw0tMVHgVWAPZM4IiKGNuGTh8TLgXcBW9nMbTqeiIheMKGTh8S6wCnA\nm2z+2XQ8ERG9YsI2mEtMAk4Hjre5qOl4IiJ6yYRNHsB/AQ8Bn2o6kIiIXjMhq60kXgC8G9gygx1G\nRAzfhCt5SCxLGdb9/Ta3Nx1PREQvmnDJA/gEcC2lvSMiIkZgQlVbSewIvAl4bu7niIgYuQlT8qi9\nq04EDrG5q+l4IiJ62YRJHsA7KXOH/LjpQCIiet2EGNtKYnXgBuClNteOfWQREd1pxGMCTpDkcQIw\nz+aQBsKKiOhaGRhxEBJbAK8Hntl0LBER48VEaPM4Fvi4zT1NBxIRMV6M65JHLXU8F9i96VgiIsaT\n8V7y+CDwFZtHmw4kImI8GbcN5hLrA1cDT7f5T7ORRUR0p8xhvrBDge8mcUREjL5xWfKQWAW4hTJq\n7t+bjisiolul5PFk7wR+kcQREdEZ47W31YHAAU0HERExXo27aivwysAdwCo285qOKSKim6XaaoEt\ngGuTOCIiOmc8Jo/nA1c2HURExHg2HpPHlsAfmw4iImI8G4/JIyWPiIgOG48N5g8Dk20eaTqeiIhu\nlwbzBW5O4oiI6KzxmDzS3hER0WHjMXmkvSMiosPGY/JIySMiosPGY4P5Kjb3NR1LREQvmBAN5pJ2\nlXSjpL9I+shA2yRxRER0Xs8kD0lLAl8FdgU2A/aV9Kxmoxo+SVObjqEdiXN0Jc7R1Qtx9kKMi6Nn\nkgewDXCz7Vm25wKn05tzk09tOoA2TW06gDZNbTqANk1tOoA2TW06gDZNbTqANkxtOoBO6qXksS5w\na8vz2+q6iIgYY72UPMZPy35ERI/rmd5WkrYDptnetT4/Ephv+9iWbXrjy0REdJGR9LbqpeQxCfgz\n8DLgdmAmsK/tGxoNLCJiAuqZaWhtz5P0XuA8YEngW0kcERHN6JmSR0REdI9eajBH0rclzZF0zSK2\n+Uq9ifBqSVuOZXwtMSwyTklvqvH9SdKlkp471jHWOIb8Pet2W0uaJ2nPsYqt3/7b+btPlXSlpGsl\nzRjD8FpjGOrvvoakcyVdVeM8YIxD7ItjfUm/knRdjeOQQbZr7FhqJ8ZuOI7a/S3rto0dR8P4m7d/\nHNnumQfwIspMgdcM8vqrgOl1eVvgsi6N84XAKnV5126Ns26zJHARcA7w+m6ME1gVuA5Yrz5fo0vj\nnAZ8ui9G4G5gUgNxrgU8ry6vSGlLfFa/bRo9ltqMsfHjqJ0462uNHkdt/p7DOo56quRh+2LgnkVs\nshtwUt32cmBVSWuORWythorT9u9s31ufXg6sNyaBLRzHUL8nwPuAHwH/6nxEA2sjzjcCP7Z9W93+\nrjEJrJ824vwnsHJdXhm42/a8jgfWj+07bF9Vlx8AbgDW6bdZo8dSOzF2w3HU5m8JDR9HbcY5rOOo\np5JHGwa6kbCRE/MwHARMbzqIgUhal3IX/4l1Vbc2kE0BJtdi+RWS9ms6oEF8A3i2pNuBq4FDG44H\nSRtRSkuX93upa46lRcTYqvHjaLA4u+04WsTvOazjqGd6Ww1D//7K3XrCQ9JLgLcCOzQdyyC+BBxh\n25LEwr9tt1iKMnf9y4Dlgd9Jusz2X5oNayFHAVfZnippY+ACSVvYvr+JYCStSLkaPrRejS60Sb/n\nY34stRFjVxxHQ8TZNcfREHEO6zgab8ljNrB+y/P16rquUxv3vgHsanuoqqOmbAWcXv6/swbwSklz\nbZ/dbFgLuRW4y/bDwMOSfgNsAXRb8tgeOAbA9l8l3QI8A7hirAORtBTwY+B7ts8aYJPGj6U2YuyK\n46iNOLviOGojzmEdR+Ot2ups4C3wxB3p/7E9p9mQFiZpA+AnwJtt39x0PIOx/XTbT7P9NMrVyru7\nMHEA/BTYUdKSkpanNPBe33BMA7kReDlAbT94BvC3sQ6iXv1+C7je9pcG2azRY6mdGLvhOGonzm44\njtr8mw/rOOqpkoek04CdgDUk3QocTSlqYfvrtqdLepWkm4EHgQO7MU7gY8BqwIn1amSu7W26MM6u\n0Mbf/UZJ5wJ/AuYD37A95smjjd/zU8B3JF1NuXD7sO1/j3WclOqdNwN/ktQ3bfNRwAZ9sXbBsTRk\njHTHcdROnN2gnb/5sI6j3CQYERHDNt6qrSIiYgwkeURExLAleURExLAleURExLAleURExLAleURE\nxLAleURUko5qcN8DDr0R0a1yn0dEJel+2yv12r4lTWpidN6Y2FLyiAlJ0pl15NBrJb1d0qeB5epE\nOKdI2lDSjZK+I+nPkr4vaec66dBNkraunzNZ0ll1UqLfSdq8rt+pftaVkv4oacU60c5vJJ1TP/vE\nOmxEX0yfVJko6neSnlrXPUXSjyTNrI/t6/ppNc5LgJNUJppaaLuIjhmryUjyyKObHsBq9d/lgGuA\nycD9La9vBMwFnk0ZBfUK4Fv1td2AM+vyccB/1+WXAFfW5bOBF9bl5SmTAU0FHq6fvQRwPnViIMpw\nEK+uy8cCH63LpwI71OUNKGMTQZlY6vfAMovaLo88OvXoqbGtIkbRoZL2qMvrUeYy6O8W29cBSLoO\nuLCuv5aSAKCMGbQngO1fSVpd0krApcD/Svo+8BPbs2shY6btWfUzTwN2pIx0+pjtn9fP/APwirr8\ncuBZLQWUlSStQBke/Wzbjy5iu+VtPzS8nyWiPUkeMeFImkqZs2A7249I+hWw7ACbPtqyPB94rGW5\n9dhZaN4L28dKOgd4NXCppF36Xuv3vvl1eW6/fU1q2WZb24/R+saSJFoTw4DbRXRK2jxiIloZuKcm\njmcB29X1cyUN94LqYuBN8ERS+pftByRtbPs625+lVC89o26/jaSNJC0B7A1cMsTnnw8c0vdE0hZt\nbve8YX6PiGFJ8oiJ6FxgkqTrKcOk/66u/z/KkNWnUEoI/bsieoDlacBWdZj1TwH71/WHSrqmrn8M\n+EVd/3vgq5R5Ev5q+8xBPrvv+SHAC2qD/HXAOweJp/927xjiN4hYLOmqGzFGasnkg7Zf23QsEYsr\nJY+IsTNQaSaiJ6XkERERw5aSR0REDFuSR0REDFuSR0REDFuSR0REDFuSR0REDFuSR0REDNv/ByZO\nRxvEva2bAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10c643650>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from pint import UnitRegistry\n",
    "ur = UnitRegistry()\n",
    "\n",
    "# ideal gas parameters for air\n",
    "gamma = 1.4\n",
    "R = 8.314 * ur.J/(ur.degK*ur.mol)\n",
    "MW = 28.966 * ur.grams/ur.mol\n",
    "\n",
    "# ambient conditions\n",
    "Pa = 1.0 * ur.atm\n",
    "Ta = ur.Quantity(20,ur.degC).to(ur.degK)\n",
    "\n",
    "def CD(P, Pa = 1.0 * ur.atm):\n",
    "    return 0.9 - 0.3*Pa/P\n",
    "\n",
    "def Qflow(T, P, Pa = 1.0 * ur.atm):\n",
    "    rho = MW*P/(R*T)\n",
    "    r = ((gamma+1)/2)**(gamma/(gamma-1))\n",
    "    if P <= r*Pa:\n",
    "        Pr = Pa/P\n",
    "        b = 2.0*(gamma/(gamma-1))*rho*P\n",
    "        Q = np.sqrt(b*(Pr**(2.0/gamma)-Pr**((1.0+1.0/gamma))))\n",
    "    else:\n",
    "        a = gamma*P*rho\n",
    "        Q = np.sqrt(a*(2.0/(gamma+1))**((gamma+1)/(gamma-1)))\n",
    "    return Q.to(ur.kg/ur.m**2/ur.sec)\n",
    "\n",
    "def Mflow(T, P, Pa, A = 1.0*ur.m**2):\n",
    "    return (CD(P,Pa)*A*Qflow(T,P,Pa)).to_base_units()\n",
    "\n",
    "P = np.linspace(1.0,2.5,100) * ur.atm\n",
    "M = [Mflow(Ta,p,Pa,A = 100*ur.cm**2) for p in P]\n",
    "\n",
    "plt.plot(P,[m.magnitude for m in M])\n",
    "plt.title('Mass Flow from an 100 cm**2 Orifice')\n",
    "plt.ylabel(M[0].units)\n",
    "plt.xlabel(P.units)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}