{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "I'm following the explanations from this blog post: [http://blog.wolfram.com/2013/10/08/a-response-to-falling-with-helium/](http://blog.wolfram.com/2013/10/08/a-response-to-falling-with-helium/).\n",
      "\n",
      "The central question to answer is: \n",
      "\n",
      ">how long it would be necessary for someone to fall in order to jump out of an airplane, fill a large balloon with helium while falling, and land safely.\n",
      "\n",
      "So roughly, my first analysis is:\n",
      "\n",
      "- jumping from an airplane: that's the free-fall problem\n",
      "- filling a large balloon with helium while falling: mmmh, here we'll likely run into trouble due to the fact that the air pressure will change as a function of altitude and thus filling a ballon with helium is not the way we usually experience it with constant air pressure at the ground level\n",
      "- land safely: this would mean that the buoyancy of the helium ballon should be a little bit less than the mass of the person falling; however buoyancy is a function of the fluid density displaced and this too varies with altitude\n",
      "\n",
      "Things that need to be determined, according to the Wolfram blog post are:\n",
      "\n",
      "- initial and final weights of the person\n",
      "- the balloon\n",
      "- the helium\n",
      "- the helium tank carried by the person falling from the plane\n",
      "\n",
      "Idea behind the calculation: we want to achieve neutral buoyancy for the falling guy. \n",
      "\n",
      "$$\n",
      "\\text{ballon buoyancy force} = \\text{weight of helium inside ballon} + \\text{weight of ballon} + \\text{weight of person}\n",
      "$$\n",
      "\n",
      "This translates to the equation used by the guy in the article:\n",
      "\n",
      "$$\n",
      "\\rho_{air} \\frac{4}{3} \\pi \\; r^3 = \\rho_{He} \\frac{4}{3} \\pi \\; r^3 + s_{balloon} 4 \\pi  \\; r^2 + m_{person}\n",
      "$$\n",
      "\n",
      "where $s_{ballon}$ is the surface density of the ballon. This makes sense because the mass of the ballon is distributed on its outer surface when fully inflated. Thus using the value of 0.00609614 grams per square centimeter calculated from a typical inflatable balloon is a good starting point.\n",
      "\n",
      "The next step in this article is to solve that equation in $r$ assuming a given air and helium density. He's saying that air density varies from 1.225 kg/m^3 at sea level to 1.1 kg/m^3 at 5000 feet (1.5 km) but we're taking it equal to 1.2 kg/m^3 as we want to be close to the surface. The density of Helium is taken constant at 1.785*10^-4 grams per cubic centimeter, or 0.18 kg/m^3 (ten times lighter than air)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Populating the interactive namespace from numpy and matplotlib\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "rho_air = 1.2 #kg/m^3\n",
      "rho_helium = 0.18 #kg/m^3\n",
      "s_ballon = 0.061 #kg/m^2\n",
      "m_person = 70 #kg"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "To find $r$ we can use the root finding method from the scipy.optimize package as discussed [in this help page](http://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html#root-finding)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.optimize import root\n",
      "def func(r):\n",
      "    return (rho_air - rho_helium) * 4./3. * pi * r ** 3 - s_ballon * 4 * pi * r ** 2 - m_person\n",
      "sol = root(func, 1)\n",
      "r = sol.x\n",
      "print \"%.2f m\" % r"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "2.60 m\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The radius of the balloon will be 2.6 meters. This represents a volume of:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "helium_volume = 4/3. * pi * r ** 3\n",
      "print \"%.2f m^3\" % helium_volume # in m^3"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "73.71 m^3\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The final mass of the jumper will thus be:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "final_mass = m_person + rho_helium * helium_volume\n",
      "print '%.2f kg' % final_mass"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "83.27 kg\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "However, the jumper in its initial state is carrying the helium with canisters that weigh 49 kg for a volume of 43 liters. How many canisters are needed to fit the helium inside? Given that each canister has a maximal design pressure, we can try to fit is much helium inside until it reaches that limit. The [HP250 model](http://www.airgas.com/content/details.aspx?id=7000000000234) has a service pressure of 156 bar. Using the ideal gas law, we know that the pressure volume product is conserved between the situations where the helium is compressed in the canisters and at ground pressure in the inflated balloon. In other words:\n",
      "$$P_{landing} V_{landing} = P_{canister} V_{canister}$$\n",
      "We take a landing pressure of 1 bar, a landing volume of 74 m^3 and a canister pressure of 156 bar. Thus the canister helium volume is:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "canister_helium_volume = 1 * helium_volume / 156.\n",
      "print \"%.2f m^3\" % (canister_helium_volume)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.47 m^3\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We therefore need, to contain the volume of helium needed to inflate the balloon, approximately"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print 0.47 / 43e-3\n",
      "number_canisters = int(0.47 / 43e-3) + 1 \n",
      "print \"%i canisters\" % (number_canisters)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10.9302325581\n",
        "11 canisters\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This leads to an initial mass of:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "initial_mass = final_mass + number_canisters * 49\n",
      "print '%.2f kg' % initial_mass"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "622.27 kg\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The next step in our problem solving is to determine the discharge time for the gas canisters. We will use the formula provided in the blog post:\n",
      "$$\n",
      "t = (F^{(1-k)/2} - 1) \\frac{2}{k - 1} \\frac{V}{C \\; A} \\left( (\\frac{k \\; P_0}{\\rho_0})(\\frac{2}{k + 1})^{(k+1)/(k-1)} \\right)^{-1/2}\n",
      "$$ \n",
      "\n",
      "> where $F$ is the fraction of the gas remaining in the container, $k$ is the heat capacity ratio (1.667 for helium), $C$ is the discharge coefficient, $A$ is the cross-sectional area of the opening, $V$ is the volume of the container, $\\rho_0$ is the initial gas density, and $P_0$ is the initial pressure. \n",
      "The discharge coefficient is the ratio of the actual discharge to the theoretical discharge. For this example, we will assume the value is 0.72."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can calculate the initial density by applying the ideal gas law:\n",
      "\n",
      "$$\n",
      "P \\; V = n \\; R \\; T\n",
      "$$\n",
      "\n",
      "Contrary to what is done in the blog post, I'd take the initial pressure to be equal to the service pressure of the canister, 156 bar. Thus the initial density would be calculated by dividing the mass of the Helium in one canister (one tenth of the final volume for which we know the density at sea level) by the volume of the canister:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "canister_helium_density = helium_volume / number_canisters * rho_helium / 43e-3\n",
      "print '%.2f kg/m^3' % canister_helium_density"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "28.05 kg/m^3\n"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This is consistent with my expectations in the sense that the density of the helium is a lot higher in the canister than in the balloon because it is compressed (which is the whole point of putting it in a canister).\n",
      "\n",
      "We can now calculate the time it takes to empty one canister (containing one tenth of the total helium volume needed) up to 90% with the formula."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "C = 0.72\n",
      "k = 1.667\n",
      "P_0 = 15.6e6 # Pa\n",
      "rho_0 = 28.23 # kg / m^3\n",
      "V = 43e-3 # m^3\n",
      "A = 3.141e-6 # m^2"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def time_to_empty(F, P_0, rho_0, V, A):\n",
      "    return (F ** ((1-k)/2.) - 1) * (2/(k-1)) * V / C / A * 1/sqrt((k * P_0 / rho_0) * (2 / (k + 1)) ** ((k + 1) / (k - 1)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "time_to_empty(0.1, P_0, rho_0, V, A) "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 12,
       "text": [
        "122.00291882678938"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fractions = arange(0, 1, 0.01)\n",
      "t = map(lambda s: time_to_empty(s, P_0, rho_0, V, A), fractions)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "-c:2: RuntimeWarning: divide by zero encountered in double_scalars\n"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot(fractions, t)\n",
      "xlabel('fraction of gas remaining in the container')\n",
      "ylabel('time needed to empty container')\n",
      "grid(True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEMCAYAAAAxoErWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtcVHX++PHXgKjlBbyBCtrgXQQZkFXSTMi7razpZktm\n4K3Wvupa7ZZWlvVrlbbMDVt33dYSzUJz17Sbq5lU3iJRshVTQ1BEUBHwSoHw+f1xYg4oOAPDXJD3\n8/E4D+YMM+e858143n4+n3M+x6CUUgghhBBVcHN2AEIIIVyXFAkhhBDVkiIhhBCiWlIkhBBCVEuK\nhBBCiGpJkRBCCFEtuxeJ0tJSQkJCGDt2LAD5+fkMHz6cHj16MGLECAoLC82vXbx4Md27d6dXr15s\n3brV3qEJIYSwwO5F4o033iAgIACDwQBAXFwcw4cP5+jRowwdOpS4uDgA0tLSWLduHWlpaWzZsoXH\nHnuMsrIye4cnhBDiJuxaJE6dOsWnn37K9OnTKb9mb/PmzcTExAAQExPDhx9+CMCmTZuIjo7Gw8MD\no9FIt27dSE5Otmd4QgghLLBrkXj88cd59dVXcXPTd3PmzBl8fHwA8PHx4cyZMwCcPn0aPz8/8+v8\n/PzIzs62Z3hCCCEsaGSvDX/88cd4e3sTEhJCUlJSla8xGAzmbqjqfm/Nc0IIISyrzSxMdmtJ7N69\nm82bN+Pv7090dDRffPEFkydPxsfHh9zcXABycnLw9vYGwNfXl6ysLPP7T506ha+vb5XbVkrJohQv\nvPCC02NwlUVyIbmQXNx8qS27FYlFixaRlZVFRkYGiYmJ3HPPPaxZs4aoqCgSEhIASEhIYNy4cQBE\nRUWRmJhIcXExGRkZHDt2jP79+9srvFtCZmams0NwGZILneRCJ7mwnd26m65X3k00b948Jk6cyMqV\nKzEajaxfvx6AgIAAJk6cSEBAAI0aNWL58uXStSSEEE5mULa0Q5zAYDDY1HS6lSQlJREREeHsMFyC\n5EInudBJLnS1PXZKkRBCiAagtsdOmZajHqvurLGGSHKhk1zoJBe2q7dF4to1eOklZ0chhBC3tnrb\n3aQUeHhAUZH2UwghRPUaXHeTwQAtW8LFi86ORAghbl31tkgAeHrChQvOjsJ5pL9VJ7nQSS50kgvb\nSZEQQghRrXo7JgEQEQELF2o/hRBCVK/BjUmAtCSEEMLe6n2RqHBjuwZH+lt1kgud5EInubBdvS8S\n0pIQQgj7qddjEs89B02awIIFTg5KCCFcnIxJCCGEqHP1ukh4eTXsIiH9rTrJhU5yoZNc2K5eF4mG\nPnAthBD2Vq/HJP77X1iyBLZudXJQQgjh4mRMQgghRJ2TIlGPSX+rTnKhk1zoJBe2kyIhhBCiWnYr\nEj/99BMDBgzAZDIREBDA/PnzAVi4cCF+fn6EhIQQEhLCZ599Zn7P4sWL6d69O7169WKrFQMNDf3s\nJrl3r05yoZNc6CQXtrPrwPXVq1e5/fbbuXbtGnfddRevvfYa27dvp0WLFjzxxBOVXpuWlsaDDz7I\nt99+S3Z2NsOGDePo0aO4uVWuYxUHX5SCxo3hyhXtpxBCiKq55MD17bffDkBxcTGlpaW0atUKoMpA\nN23aRHR0NB4eHhiNRrp160ZycvJNt28wNOwuJ+lv1UkudJILneTCdo3sufGysjJCQ0NJT09n5syZ\n9OnThw0bNrBs2TJWr15NWFgYS5YswcvLi9OnTxMeHm5+r5+fH9nZ2VVuNzY2FqPRCIBSXmzdamLS\npAhA/1KUNzNlvWGsl3OVeJy5npqa6lLxOHM9NTXVpeJx5HpSUhKrVq0CMB8va8Mh10lcuHCBkSNH\nEhcXR0BAAO3atQNgwYIF5OTksHLlSmbPnk14eDiTJk0CYPr06YwZM4bx48dXDvi6JlNoKLz1FvTr\nZ+9PIYQQ9ZdLdjeV8/T05N5772Xfvn14e3tjMBgwGAxMnz7d3KXk6+tLVlaW+T2nTp3C19fXim03\n3O4mIYSwN7sViby8PAp/mTOjqKiIbdu2ERISQm5urvk1GzduJCgoCICoqCgSExMpLi4mIyODY8eO\n0b9/f4v78fJquFNzXN/V0pBJLnSSC53kwnZ2G5PIyckhJiaGsrIyysrKmDx5MkOHDuXhhx8mNTUV\ng8GAv78/K1asACAgIICJEycSEBBAo0aNWL58OQaDweJ+pCUhhBD2c9MxidLSUuLj43n88ccdGdNN\nXd+v9oc/gL8/zJ3rxKCEEMLF2WVMwt3dnffee6/WQTmCtCSEEMJ+LI5J3HXXXcyaNYuvv/6a/fv3\nmxdX0ZCLhPS36iQXOsmFTnJhO4tjEgcOHMBgMPD8889Xen7Hjh12C6omPD0hLc3ZUQghxK2pXt9P\nAmDDBnj/ffj3v50YlBBCuDi7XSeRm5vLtGnTGDVqFKDNsbRy5cqaR2gnDbm7SQgh7M1ikYiNjWXE\niBGcPn0agO7du7N06VK7B2athlwkpL9VJ7nQSS50kgvbWSwSeXl5PPDAA7i7uwPg4eFBo0Z2nfKp\nRhpykRBCCHuzWCSaN2/O+fPnzet79+7F09PTrkHVREMuEuWTegnJRUWSC53kwnYWmwRLlixh7Nix\nHD9+nIEDB3Lu3Dk2bNjgiNis0pCn5RBCCHuz6uymkpISjhw5AkDPnj3x8PCwe2DVqWqEvkkTrTXR\ntKmTgnKSpKQk+Z/SLyQXOsmFTnKhq+3ZTVYNLiQnJ5OZmcm1a9fMF9I9/PDDNd6ZvZR3OTW0IiGE\nEPZmsSXx0EMPcfz4cUwmk3nwGmDZsmV2D64qVVXD7t3hk0+gRw+nhCSEEC7Pbi2JlJQU0tLSrJqR\n1Vka8uC1EELYk8WzmwIDA8nJyXFELLXWUIuEnAOuk1zoJBc6yYXtLLYkzp07R0BAAP3796dJkyaA\n1mzZvHmz3YOzlpzhJIQQ9mFxTKK6SuysMwaq6lebOhUGDYJp05wSkhBCuDy7jUnUh9PHGmp3kxBC\n2Fu1YxKDBg0CtCuuW7RoUWlp2bKlwwK0RkMtEtLfqpNc6CQXOsmF7aptSezatQuAy5cvOyyY2vL0\nhBMnnB2FEELceiye3VTu7NmznDx50rxY8tNPPzFgwABMJhMBAQHMnz8fgPz8fIYPH06PHj0YMWIE\nhRVGnBcvXkz37t3p1asXW7dutfpDeHo2zIHr+tAV6CiSC53kQie5sJ3FIrF582a6d++Ov78/Q4YM\nwWg0Mnr0aIsbbtq0KTt27CA1NZWDBw+yY8cOdu7cSVxcHMOHD+fo0aMMHTqUuLg4QLtPxbp160hL\nS2PLli089thjlJWVWfUhvLwaZneTEELYm8Ui8dxzz7Fnzx569OhBRkYG27dvZ8CAAVZt/Pbbbweg\nuLiY0tJSWrVqxebNm4mJiQEgJiaGDz/8EIBNmzYRHR2Nh4cHRqORbt26kZycbNV+ZExCSC50kgud\n5MJ2Fs9u8vDwoG3btpSVlVFaWkpkZCR/+MMfrNp4WVkZoaGhpKenM3PmTPr06cOZM2fw8fEBwMfH\nhzNnzgBw+vRpwsPDze/18/MjOzu7yu3GxsZiNBoB8PLy4rbbTFy4EAHoX4ryZqasN4z1cq4SjzPX\nU1NTXSoeZ66npqa6VDyOXE9KSmLVqlUA5uNlbVi8TmLYsGFs3LiR+fPnk5eXh7e3N/v27WP37t1W\n7+TChQuMHDmSxYsXM378eAoKCsy/a926Nfn5+cyePZvw8HAmTZoEwPTp0xkzZgzjx4+vHHAV5/oe\nOwajR8OPP1odkhBCNCh2u8f1pk2buP3221m6dCmjRo2iW7dufPTRRzXaiaenJ/feey8pKSn4+PiQ\nm5sLQE5ODt7e3gD4+vqSlZVlfs+pU6fw9fW1cvsNs7tJCCHszWKReOmll3B3d8fDw4PY2FjmzJnD\nX/7yF4sbzsvLM5+5VFRUxLZt2wgJCSEqKoqEhAQAEhISGDduHABRUVEkJiZSXFxMRkYGx44do3//\n/lZ9iPKzm2pRJOu167taGjLJhU5yoZNc2M5ikajqVNRPP/3U4oZzcnK45557MJlMDBgwgLFjxzJ0\n6FDmzZvHtm3b6NGjB1988QXz5s0DICAggIkTJxIQEMDo0aNZvny51TPPNmkCjRpBUZFVLxdCCGGl\nasck/v73v7N8+XLS09Pp2rWr+flLly4xaNAg1q5d67AgK6quX619ezhwADp0cEJQQgjh4mo7JlFt\nkbhw4QIFBQXMmzePV155xbzxFi1a0KZNG9uitUF1H7RnT9i0CXr1ckJQQgjh4up84NrT0xOj0Uhi\nYiJ+fn40btwYNzc3rly5YtUV147WEAevpb9VJ7nQSS50kgvbWbxOYtmyZbz44ot4e3tXun3p999/\nb9fAaqohFgkhhLA3i9dJdO3aleTkZKd2MVVUXZNp4kS47z6IjnZCUEII4eLsdp1E586dXW5q8Kp0\n6gQVLrMQQghRByx2N/n7+xMZGcm9995L48aNAa0iPfHEE3YPribuuAOOHHF2FI6VlJRkvhy/oZNc\n6CQXOsmF7SwWic6dO9O5c2eKi4spLi5GKWX19QuOdMcdUIPZxYUQQljB4phEuUuXLgHaKbDOVF2/\n2nffwUMPgYuNpwshhEuw25jE999/T0hICH369KFPnz7069eP//3vf7UK0p7uuEO7O11Dm5pDCCHs\nyWKReOSRR3j99dfNd6RbsmQJjzzyiCNiqxEvLzAYoMIEs7c8OQdcJ7nQSS50kgvbWSwSV69eJTIy\n0rweERHBlStX7BpUbZW3JoQQQtQNi2MS48aNo1+/fkyePBmlFGvXriUlJYWNGzc6KsZKbtavNnYs\nTJsGv0wsK4QQ4hd2G5N4++23OXv2LOPHj2fChAmcO3eOt99+u1ZB2pvRKC0JIYSoSxaLROvWrVm2\nbBn79+9n//79vPHGG7Rq1coRsdVYQ+tukv5WneRCJ7nQSS5sZ7FIDBs2zHzzIID8/HxGjhxp16Bq\n6447IDPT2VEIIcStw+KYhMlkMt9M/GbPOcrN+tWSk2HmTEhJcXBQQgjh4uw2JuHu7s6JCn04mZmZ\nuLlZfJtTNLTuJiGEsDeLR/s///nPDB48mMmTJ/PQQw9x9913s2jRIkfEVmPe3nD1KrjoGbp1Tvpb\ndZILneRCJ7mwncW5m0aNGkVKSgp79+7FYDCwdOlS2rVr54jYasxggM6dtdZEQICzoxFCiPrPqn6j\ndu3aMXbsWH7961/XqEBkZWURGRlJnz59CAwMJD4+HoCFCxfi5+dHSEgIISEhfPbZZ+b3LF68mO7d\nu9OrVy+21mLGvoY0eC2zW+okFzrJhU5yYTuLLQlbeHh4sHTpUkwmE5cvX6Zfv34MHz7cPNX49dON\np6WlsW7dOtLS0sjOzmbYsGEcPXq0RmMgMi4hhBB1x64j0O3bt8dkMgHQvHlzevfuTXZ2NkCVo+yb\nNm0iOjoaDw8PjEYj3bp1Izk5uUb7bEhFQvpbdZILneRCJ7mwncWWxBNPPMG0adPo06ePTTvKzMzk\nwIEDhIeHs2vXLpYtW8bq1asJCwtjyZIleHl5cfr0acLDw83v8fPzMxeVimJjYzEajQB4eXlhMpnM\nzcqioiS+/RZAWy//kpT/XtZvzfVyrhKPM9dTU1NdKh5nrpefqu8q8ThyPSkpiVWrVgGYj5e1YfE6\nibfeeotVq1ZRUlLC1KlTiY6OxtPTs0Y7uXz5MhERETz33HOMGzeOs2fPmsc2FixYQE5ODitXrmT2\n7NmEh4czadIkAKZPn86YMWMYP368HrCFc3137oSnnoLdu2sUohBC3NLsdp3EjBkz2LVrF6tXryYz\nM5OgoCAefPBBduzYYdUOSkpKmDBhAg899BDjfpl5z9vbG4PBgMFgYPr06eYuJV9fX7Iq3Kj61KlT\n+Pr61ugDNaSBayGEsDerxiRKS0v54YcfOHz4MO3atSM4OJjXX3+dBx544KbvU0oxbdo0AgICmDt3\nrvn5nJwc8+ONGzcSFBQEQFRUFImJiRQXF5ORkcGxY8fo379/jT5Qx45w/jz8/HON3lYvXd/V0pBJ\nLnSSC53kwnYWxyQef/xxPvroI+655x6effZZ80H76aefpmfPnjd9765du3j33Xfp27cvISEhACxa\ntIj333+f1NRUDAYD/v7+rFixAoCAgAAmTpxIQEAAjRo1Yvny5TW+n7a7u1YosrKgW7cavVUIIcR1\nLI5JvPPOO0ycOJFmzZrd8LvCwkK8vLzsFlxVrOlXi4iABQtg6FDHxCSEEK7ObmMSa9asuaFADP3l\n6OvoAmEtua+EEELUjWqLRFFREefPnycvL4/8/HzzkpmZWeVpqa6koQxeS3+rTnKhk1zoJBe2q3ZM\nYsWKFbzxxhucPn2afv36mZ9v0aIFs2bNckhwtXXHHfDll86OQggh6j+LYxLx8fHMmTPHUfFYZE2/\n2u7d8Ic/8MtFdUIIIWo7JmGxSBQVFbF8+XJ27tyJwWBg8ODBzJw5k6ZNm9Y6WFtY80EvXwYfHygs\nBA8PBwUmhBAuzG4D1w8//DBpaWnMmTOHWbNmcejQISZPnlyrIB2leXPo1AkOH3Z2JPYl/a06yYVO\ncqGTXNjO4nUShw4dIi0tzbx+zz33EFAPbtYQGgoHDkDfvs6ORAgh6i+LLYnQ0FD27NljXt+7d2+l\ngWxXFRKiFYlbWfmkXkJyUZHkQie5sJ3FlsS+ffsYNGgQnTp1wmAwcPLkSXr27ElQUBAGg4GDBw86\nIs4aCwmBjz5ydhRCCFG/WRy4zrRwwYEtU9DWhrWDL+fPQ5cuUFAAbna9a4bzJCUlyf+UfiG50Eku\ndJILXW0Hri22JIxGIwUFBWRlZXHt2jXz86GhoTXemSO1aQNeXnD8uMzhJIQQtWWxJbFgwQJWrVpF\nly5dKt1G1NqpwutaTarhffdBdDRMnGjnoIQQwsXZrSWxbt060tPTady4ca0Cc6bywWspEkIIUTsW\ne+v79OlDQUGBI2Kpc7f6GU5yDrhOcqGTXOgkF7az2JJ45plnCAkJITAwkCZNmgBas2Xz5s12D85W\noaGwfz8oBTW8LYUQQgisGJPo3bs3M2fOJDAw0DwmYTAYGDJkiEMCvF5N+tWU0qbnOHAAangXVCGE\nuKXYbUyiefPmLjXBX00YDPqV11IkhBCi5iyOSQwePJj58+ezZ88e9u/fb17qi5AQrcvpViT9rTrJ\nhU5yoZNc2M5iS2L//v0YDAb27t1b6XlnnQJbUyEh8P77zo5CCCHqJ4tjErbIysri4Ycf5uzZsxgM\nBh555BHmzJlDfn4+DzzwACdOnMBoNLJ+/XrzrVAXL17M22+/jbu7O/Hx8YwYMaJywDXsV/vxR+1e\n13I7UyFEQ2a3qcJzc3OZNm0ao0aNAiAtLY2VK1datXEPDw+WLl3KoUOH2Lt3L3/72984fPgwcXFx\nDB8+nKNHjzJ06FDi4uLM2163bh1paWls2bKFxx57jLKyshp/qIq6dIFLlyAnx6bNCCFEg2SxSMTG\nxjJixAhOnz4NQPfu3Vm6dKlVG2/fvj0mkwnQBsB79+5NdnY2mzdvJiYmBoCYmBg+/PBDADZt2kR0\ndDQeHh4YjUa6detGcnJyrT5YOTc3iIyEzz+3aTMuSfpbdZILneRCJ7mwncUxiby8PB544AHz//Y9\nPDxo1Mji226QmZnJgQMHGDBgAGfOnMHHxwcAHx8fzpw5A8Dp06cJDw83v8fPz4/s7OwbthUbG2ue\nWNDLywuTyWSexKv8S1Fx3WiEbdsimDy56t/Lev1fL+cq8ThzPTU11aXiceZ6amqqS8XjyPWkpCRW\nrVoF2DgRq7JgyJAhKi8vT5lMJqWUUnv27FF33323pbdVcunSJRUaGqo2btyolFLKy8ur0u9btWql\nlFJq1qxZ6t133zU/P23aNPXvf/+70mutCPkG6elKtW+vVFlZjd8qhBC3hNocO5VSymKTYMmSJYwd\nO5bjx48zcOBAzp07x4YNG6wuQiUlJUyYMIHJkyczbtw4QGs95Obm0r59e3JycvD29gbA19eXrKws\n83tPnTqFbx1c4NClC9x+O/zvfxAUZPPmhBCiwbA4JtGvXz++/PJLdu3axYoVKzh06BDBwcFWbVwp\nxbRp0wgICGDu3Lnm56OiokhISAAgISHBXDyioqJITEykuLiYjIwMjh07Rv/+/WvzuW4wYgRs21Yn\nm3IZ13e1NGSSC53kQie5sJ1VgwseHh4EBgbWeOO7du3i3XffpW/fvoSEhADaKa7z5s1j4sSJrFy5\n0nwKLEBAQAATJ04kICCARo0asXz5cgx1NOnS8OHwz3/CE0/UyeaEEKJBsOt1EvZQ23N9CwuhUyc4\ndw6aNrVDYEII4cLsdp3ErcLLC/r0gd27nR2JEELUH1YViU2bNvHkk0/y5JNP8tFHH9k7JrsZMQK2\nbnV2FHVH+lt1kgud5EInubCdxSIxb9484uPj6dOnDwEBAcTHxzN//nxHxFbnhg+/9QavhRDCniyO\nSQQFBZGamoq7uzsApaWlmEwmvv/+e4cEeL3a9qsBlJRA27bafE7t2tVxYEII4cLsNiZhMBgoLCw0\nrxcWFtbZGUeO5uEBQ4bcWl1OQghhTxaLxPz58wkNDSU2NpaYmBj69evHM88844jY7GLCBFi3ztlR\n1A3pb9VJLnSSC53kwnYWr5OIjo5myJAhfPvttxgMBuLi4ujQoYMjYrOL++6DOXPg/Hlo08bZ0Qgh\nhGuzOCYxdOhQtm/fbvE5R7FlTKLcAw/APffAo4/WUVBCCOHi6nxMoqioiPPnz3Pu3Dny8/PNS2Zm\nZpUzs9YnkybB2rXOjkIIIVxftUVixYoVhIWFceTIEfr162deoqKimDVrliNjrHOjRkFaWv2/W530\nt+okFzrJhU5yYbtqi8TcuXPJyMjg1VdfJSMjw7wcPHiw3heJxo21AWy597UQQtxcg5m76XpffQX/\n93/gpMs9hBDCoWTuphq66y64cEGKhBBC3EyDLRJubvDgg/V7AFv6W3WSC53kQie5sF2110mkpKSY\nmydVXWEdGhpq18AcISZGOxV24UKZPlwIIapS7ZhEREQEBoOBoqIiUlJS6Nu3LwAHDx4kLCyMPXv2\nODTQcnU1JlFu1CjtuokpU+psk0II4XLqfEwiKSmJHTt20LFjR/bv309KSgopKSkcOHCAjh072hSs\nK3niCVi6FOrX8L0QQjiGxTGJH374gaCgIPN6YGAghw8ftmtQjjR8OJSVgZMuILeJ9LfqJBc6yYVO\ncmE7i3M39e3bl+nTp/PQQw+hlOK9994jODjYEbE5hMEAjz+utSaGDXN2NEII4VosXidRVFTE3//+\nd77++msA7r77bmbOnElTK0Z6p06dyieffIK3t7f5/hMLFy7kX//6F+1+uaHDokWLGD16NACLFy/m\n7bffxt3dnfj4eEaMGHFjwHU8JgHw009gNMKOHdC7d51uWgghXEJtj51WXUx39epVTp48Sa9evWq0\n8a+//prmzZvz8MMPm4vEiy++SIsWLXjiiScqvTYtLY0HH3yQb7/9luzsbIYNG8bRo0dxc6vcI2aP\nIgHaGU45ObBiRZ1vWgghnM5uF9Nt3ryZkJAQRo0aBcCBAweIioqyauODBw+mVatWNzxfVaCbNm0i\nOjoaDw8PjEYj3bp1Izk52ar91IXHHoP16yE312G7tJn0t+okFzrJhU5yYTuLYxILFy7km2++ITIy\nEoCQkBCOHz9u006XLVvG6tWrCQsLY8mSJXh5eXH69GnCw8PNr/Hz86t2ttnY2FiMRiMAXl5emEwm\nIiIiAP1LUZv1KVPg979PYu7c2r1f1p23Xs5V4nHmempqqkvF48z11NRUl4rHketJSUmsWrUKwHy8\nrA2L3U0DBgzgm2++ISQkhAMHDgDaYPbBgwet2kFmZiZjx441dzedPXvWPB6xYMECcnJyWLlyJbNn\nzyY8PJxJkyYBMH36dMaMGcP48eMrB2yn7iaAvDzo1Qv27oVu3eyyCyGEcAq7dTf16dOHtWvXcu3a\nNY4dO8bs2bMZOHBgrYIE8Pb2xmAwYDAYmD59urlLydfXl6ysLPPrTp06ha+vb633Uxtt22pnOj33\nnEN3K4QQLstikVi2bBmHDh2iSZMmREdH07JlS/7617/Weoc5OTnmxxs3bjRfgxEVFUViYiLFxcVk\nZGRw7Ngx+vfvX+v91NbcudoMsSkpDt91jV3f1dKQSS50kgud5MJ2FsckmjVrxqJFi1i0aFGNNx4d\nHc2XX35JXl4enTp14sUXXzT3mRoMBvz9/Vnxy+lEAQEBTJw4kYCAABo1asTy5curnDPK3po1g+ef\nh3nzYNs2h+9eCCFcSrVjEmPHjtVfdF1flsFgYPPmzfaPrgr2HJMoV1ICffrAm29CFZdqCCFEvVPb\nY2e1LYknn3wS0LqEcnNzzVdcv//++/j4+NQ+0nrAwwOWLNFuSnTwINx2m7MjEkII57B4dlO/fv1I\nua6DvqrnHMURLYly998P3btDLXraHCIpKcl86ltDJ7nQSS50kgud3c5uunr1Kunp6eb148ePc/Xq\n1RrvqD6Kj4d//UtrTQghRENksSWxZcsWHnnkEfz9/QHtuod//vOfjBw50iEBXs+RLQmAt97SCsXu\n3eDu7rDdCiFEnbLr3E0//fQTR44cAaBXr140adKk5hHWEUcXibIyiIiA3/4W5sxx2G6FEKJO2a27\n6cqVK7z66qu8+eabBAcHc/LkST7++ONaBVkfublprYn/9//g0CFnR1OZnAOuk1zoJBc6yYXtLBaJ\nKVOm0LhxY3bv3g1Ax44defbZZ+0emCvp2RNeeUW7zWkDGY4RQgigBmc3VZy7KTg4mO+++84hAV7P\n0d1N5ZSChx7SLrb75z8dvnshhLCJ3bqbmjRpQlFRkXk9PT3dqWMSzmIwwD/+od2YaN06Z0cjhBCO\nYbFILFy4kFGjRnHq1CkefPBB7rnnHl555RVHxOZyWrTQCsTs2fDLOL5TSX+rTnKhk1zoJBe2szh3\n04gRIwgNDeWbb75BKUV8fDxt27Z1RGwuKTQU4uLg17/WphRv08bZEQkhhP1YHJMoKytj7dq1ZGRk\n8Pzzz3MIz8mAAAAex0lEQVTy5Elyc3OdMkMrOG9M4npPPQXJybB1KzRu7OxohBDi5ux2ncTvf/97\n3Nzc2LFjB4cPHyY/P58RI0awb9++WgdrC1cpEqWlMH68dg+Kf/1LG7MQQghXZbeB62+++Ybly5fT\ntGlTAFq3bk1JSUnNI7zFuLvD2rXafScWL3ZODNLfqpNc6CQXOsmF7SyOSTRu3JjS0lLz+rlz53Bz\ns1hbGoTmzeHTT2HwYG1Qe/ZsZ0ckhBB1y2J307vvvsv69etJSUkhJiaGDRs28PLLLzNx4kRHxViJ\nq3Q3VZSZCXffDS++CFOmODsaIYS4kV3nbjp8+DDbt28HYOjQofTu3bvmEdYRVywSoJ0SGxkJS5dq\nV2YLIYQrsduYBECPHj247777GDt2LM2aNePkyZM13tGtrmdP2LJFu0f26tWO2af0t+okFzrJhU5y\nYTuLYxLLli3jxRdfxNvbG/cKc2V///33dg2sPurbF774Qrvl6aVL2p3thBCiPrPY3dS1a1eSk5Np\nU4urxqZOnconn3yCt7e3uajk5+fzwAMPcOLECYxGI+vXr8fLywuAxYsX8/bbb+Pu7k58fDwjqrjB\ntKt2N1WUkQHDhsHUqfDMM3J6rBDC+ezW3dS5c2datmxZq6CmTJnCli1bKj0XFxfH8OHDOXr0KEOH\nDiUuLg6AtLQ01q1bR1paGlu2bOGxxx6jrKysVvt1Nn9/2LkTEhNh1iy4ds3ZEQkhRO1YLBL+/v5E\nRkayePFilixZwpIlS3j99det2vjgwYNp1apVpec2b95MTEwMADExMXz44YcAbNq0iejoaDw8PDAa\njXTr1o3k5OSafh6X0aGDVijS07UpPC5cqPt9SH+rTnKhk1zoJBe2szgm0blzZzp37kxxcTHFxcU2\n7/DMmTP4+PgA4OPjw5kzZwA4ffo04eHh5tf5+fmRnZ1d5TZiY2MxGo0AeHl5YTKZzDc7L/9SuMK6\npyf86U9JvPkmDBwYwccfw4kTrhPfrbRezlXiceZ6amqqS8XjzPXU1FSXiseR60lJSaxatQrAfLys\nDatOgbVFZmYmY8eONY9JtGrVioKCAvPvW7duTX5+PrNnzyY8PJxJkyYBMH36dMaMGcP48eMrB1wP\nxiSqsmwZvPwyvPMOjBnj7GiEEA2NXU+BrUs+Pj7k5uYCkJOTg7e3NwC+vr5kZWWZX3fq1Cl8fX0d\nHZ7dzJ4N//43PPIIPP+8NveTEEK4OocXiaioKBISEgBISEhg3Lhx5ucTExMpLi4mIyODY8eOOW2m\nWXu56y5trqedO2HkSDh92rbtXd/V0pBJLnSSC53kwnZ2LRLR0dEMHDiQI0eO0KlTJ9555x3mzZvH\ntm3b6NGjB1988QXz5s0DICAggIkTJxIQEMDo0aNZvnw5hlvw3FEfH2168cGDtXtTbNzo7IiEEKJ6\nFsckjhw5wmOPPUZubi6HDh3i4MGDbN68meeee85RMVZSX8ckqrJnD0yeDBER2nQeLVo4OyIhxK3K\nbmMSM2bMYNGiRTRurN1ZJygoiPfff7/mEYob3HknHDigXWwXGAiffebsiIQQojKLReLq1asMGDDA\nvG4wGPDw8LBrUA1Jixbw1luwcqU2jcfDD8P589a9V/pbdZILneRCJ7mwncUi0a5dO3788Ufz+oYN\nG+jQoYNdg2qIhg2D77/X7pkdEKAVjnp6wbkQ4hZicUwiPT2dRx55hN27d9OqVSv8/f1Zu3atTRdn\n2OJWGpOoTmqq1qooKYG//Q1+9StnRySEqO/sej8JgCtXrlBWVkYLJ4+uNoQiAVorYs0amD9fa2X8\n+c/QqZOzoxJC1Fd2G7guKCjgjTfe4LnnnuOZZ55h9uzZzJkzp1ZBCuu5uUFMjHYzo86dwWSC556D\nixf110h/q05yoZNc6CQXtrNYJMaMGcOJEyfo27cvYWFh9OvXj379+jkiNoE2sP3yy1oXVFYWdOsG\nf/kLXL3q7MiEEA2Bxe6m0NBQ9u/f76h4LGoo3U3VSUvTpvXYvRvmzYMZM+C225wdlRDC1dltTOK1\n116jZcuWjB07liZNmpifb926dc2jrAMNvUiU278fXnoJvvkGnnwSfv97aN7c2VEJIVyV3cYkmjZt\nyp/+9CfCw8PNXU1hYWG1ClLUndBQmDs3if/+F779VrvR0YIFcPassyNzDul71kkudJIL21ksEkuW\nLCE9PZ0TJ06QkZFBRkYGx48fd0Rswgp9+8K6dbBrF5w7Bz17wqOPwuHDzo5MCHErsNjdNGLECDZu\n3EizZs0cFdNNSXfTzZ09q11bsWIFBAdrU5SPGaOdLSWEaLjsNiYxbtw4Dh06RGRkpHlMwmAwEB8f\nX7tIbSRFwjo//6y1MOLjIT9fa11MmQK/3L5DCNHA2G1MYty4cTz77LMMHDjQPCYhp8C6hpv1tzZp\nos0D9e23kJioXW/Rsyf87nfw+ee33pQf0vesk1zoJBe2s3iP69jYWAeEIezFYID+/bXl9dfh3Xfh\nqae0SQRjY7UL9rp0cXaUQghXVW130/33388HH3xAUFDQjW8yGDh48KDdg6uKdDfVjQMHtPttJyZC\njx7afS0mToRWrZwdmRDCHup8TOL06dN07NiREydO3LBhg8HAHXfcUbtIbSRFom6VlMCWLdo8Uf/9\nLwwZonVJRUXJdRdC3ErqfEyiY8eOACxfvhyj0VhpWb58ee0jFXWmLvpbPTxg7FhYv16b9mPiRHjv\nPfD1hQkT4P33K88X5aqk71knudBJLmxnceB669atNzz36aef2rxjo9FI3759CQkJoX///gDk5+cz\nfPhwevTowYgRIygsLLR5P8J6LVvCQw/Bxx9DRoZWPNauBT8/7fHKlQ33Yj0hGqpqu5v+/ve/s3z5\nctLT0+natav5+UuXLjFo0CDWrl1r0479/f1JSUmpNL3HU089Rdu2bXnqqad45ZVXKCgoIC4urnLA\n0t3kcIWF2q1VN27UuqSCgrSiMXYs9O6tDY4LIVxbnY9JXLhwgYKCAubNm8crr7xi3niLFi1o06aN\nbdGiFYl9+/ZV2lavXr348ssv8fHxITc3l4iICH744YfKAUuRcKqffoIdO+Cjj7QWh4eHdrHe6NEQ\nEQG33+7sCIUQVbH7TYfqWpcuXfD09MTd3Z1HH32UGTNm0KpVKwoKCgBQStG6dWvzujlgKRJmSUlJ\nREREOG3/SsHBg1or47PPtEkHBw6EESO0JTDQca0MZ+fClUgudJILXW2PnRavk7CXXbt20aFDB86d\nO8fw4cPp1atXpd8bDAYM1RxhYmNjzbdP9fLywmQymb8I5QNVsm7/dYMBCgqSCA+HefMiuHAB4uOT\n+Oor+NvfIigqgsDAJEJD4f/+L4LOne0XTzlXyo+z1lNTU10qHmeup6amulQ8jlxPSkpi1apVADbd\nbtppLYmKXnzxRZo3b85bb71FUlIS7du3Jycnh8jISOluqqeUguPHYft2bfniC/D0hMhIrVsqIkI7\ng0oI4Rh2m5bDHq5evcqlS5cA7d7ZW7duJSgoiKioKBISEgBISEhg3LhxzghP1AGDAbp2hUce0eaQ\nOnNGG/gOCoING7TJB7t21eaTeucdOHZMKyxCCNfilJZERkYG9913HwDXrl1j0qRJzJ8/n/z8fCZO\nnMjJkycxGo2sX78eLy+vygFLS8IsqR73t5aVadOZf/WVtuzcCcXFcNdd2rjGwIHaPTMq3Ofqpupz\nLuqa5EInudDVqzEJf39/c19hRa1bt+bzzz93QkTC0dzcoE8fbZk5U3vu5EmtWOzZo13Q98MP2v0y\nBgyA8HDtp9Eop9wK4UguMSZRE9KSaDguX4Z9+7RbtO7dq/0sKdEnLAwL0xYfH2dHKoTrq3enwNaW\nFImGSynIztamP09O1gpISop2bUZYmNY9Vb506CAtDiEqkiLRAEl/q1Y4MjJg9eokfv45gpQU7XoN\nd3cICQGTSVuCg6F7d2jktJO+HUe+FzrJha5ejUkIUVcMBu1+GOWn1YLe4jhwQFs++ACeew5ycrRp\nRPr21ZagIG2Ru/UJUT1pSYgG49Il+P57fTl4EP73P611ERSkDaIHBOgD6hWmFROi3pPuJiFqQSmt\nhfH993DokL6kpWljHQEBWuujd2/o1Uv72bGjjHeI+keKRAMk/a26us5FeZfV4cNa0ThyRHt8+DAU\nFWl38+vZU1969NCWZs3qLIRak++FTnKhkzEJIeqQwaDdR8PPD4YPr/y7ggKtaBw5ol3LsX49HD0K\nP/6odVH16KENkpcv3bppV5ffdptzPosQtpCWhBB1pKxMu7vfsWP68uOP2s+MDGjbVisW5UWjSxd9\nadNGurCEfUl3kxAurLRU67768UdIT6+8HD+uFZguXcDfX1sqPr7jDtfoxhL1mxSJBkj6W3X1PRcF\nBVqxOH5ca3WUL5mZcOIENG+uTUlyxx3aUv64c2ftp5eX3hKp77moS5ILnYxJCFGPtWoF/fppy/WU\n0mbRPXFCKxqZmdoA+n//qz134oT2us6dtaVRI9i1Czp10hc/P2ja1JGfSNwqpCUhxC3gwgW9YGRl\n6cvJk3DqlNbV1bKlXjCuX3x9taV5c2d/EmEv0t0khKhWWRmcPasVjPIlK0srHtnZeiHx8NALRseO\nVS/t21s/hbtwHVIkGiDpb9VJLnS1zYVSUFioFYvTp7WlvIjk5OjrZ89qrZIOHW5c2rfXf7ZvDy1a\nOPesLfle6GRMQghhE4NBGxtp1QoCA6t/XVkZ5OVpRSM3VysgOTnamVs7d1Z+Tim9YPj46D+rWpo3\nl9OAXZG0JIQQdnP5slY0zpy58WfF5exZrfh4e1e9tGt340/p8qoZ6W4SQtRrly/DuXNawTh7Vise\nFdcrPs7L065gb9dOX9q2rfyz/HH50rJlw26pSJFogKS/VSe50DWEXCilndF17tyNS16etpw7B8eP\nJ1FcHEFeHvz0k3Zle5s2euEoXy9/ruJ6mzba9Sfu7s7+tHXjlhmT2LJlC3PnzqW0tJTp06fz9NNP\nOzskl5WamnrLHwysJbnQNYRcGAzaAdzLS5sfqzp//Wsqc+dGAPDzz3D+vLacO6c/zsvTzu767jv9\nufLl4kWtBdKmjTYvV/nP6pZWrfSft8oNrlzqY5SWljJr1iw+//xzfH19+dWvfkVUVBS9e/d2dmgu\nqbCw0NkhuAzJhU5yoauYiyZN9NN4rVVaqp3xlZ+vF46CAv3xsWPa7/LztefLHxcWalPNlxcMa5fy\nwufhYYdk1JJLFYnk5GS6deuG0WgE4He/+x2bNm2SIiGEcAp3d73r6WYtluuVlWmtkIKC6peTJ298\nrrBQW5o21YpFxcLh5QWenpXXKz5X8WddDuq7VJHIzs6mU6dO5nU/Pz+++eYbJ0bk2jIzM50dgsuQ\nXOgkFzpn5cLNTT+I+/vX7L1KaYP45QWjsFArIBcu6I9zcrSpWcqfL/9d+ePy/X//vTaAbwuXKhIG\nK089sPZ1DUFCQoKzQ3AZkgud5ELXUHNx5kzd3L/dpYqEr68vWVlZ5vWsrCz8/PwqvUbObBJCCMdx\nc3YAFYWFhXHs2DEyMzMpLi5m3bp1REVFOTssIYRosFyqJdGoUSPefPNNRo4cSWlpKdOmTZNBayGE\ncCKXakkAjB49miNHjvDmm2+SkJBA9+7deeWVV6p87Zw5c+jevTvBwcEcOHDAwZE6zpYtW+jVq1e1\nuVi7di3BwcH07duXQYMGcfDgQSdE6RiWclHu22+/pVGjRvznP/9xYHSOZU0ukpKSCAkJITAw8Ja+\ndsJSLvLy8hg1ahQmk4nAwEBWrVrl+CAdYOrUqfj4+BAUFFTta2p83FQu6Nq1a6pr164qIyNDFRcX\nq+DgYJWWllbpNZ988okaPXq0UkqpvXv3qgEDBjgjVLuzJhe7d+9WhYWFSimlPvvsswadi/LXRUZG\nqnvvvVdt2LDBCZHanzW5KCgoUAEBASorK0sppdS5c+ecEardWZOLF154Qc2bN08ppeWhdevWqqSk\nxBnh2tVXX32l9u/frwIDA6v8fW2Omy7XkoDK10t4eHiYr5eoaPPmzcTExAAwYMAACgsLOXPmjDPC\ntStrcnHnnXfi6ekJaLk4deqUM0K1O2tyAbBs2TJ++9vf0s7Wc/9cmDW5eO+995gwYYL55I+2bds6\nI1S7syYXHTp04OLFiwBcvHiRNm3a0OhWuSS6gsGDB9OqVatqf1+b46ZLFomqrpfIzs62+Jpb8eBo\nTS4qWrlyJWPGjHFEaA5n7fdi06ZNzJw5E7h1T5e2JhfHjh0jPz+fyMhIwsLCWLNmjaPDdAhrcjFj\nxgwOHTpEx44dCQ4O5o033nB0mC6hNsdNlyyl1v7DVtedDnsrHhBq8pl27NjB22+/za5du+wYkfNY\nk4u5c+cSFxdnnszs+u/IrcKaXJSUlLB//362b9/O1atXufPOOwkPD6d7TS4drgesycWiRYswmUwk\nJSWRnp7O8OHD+e6772jRooUDInQtNT1uumSRsOZ6ietfc+rUKXx9fR0Wo6NYkwuAgwcPMmPGDLZs\n2XLT5mZ9Zk0uUlJS+N3vfgdog5WfffYZHh4et9yp1NbkolOnTrRt25bbbruN2267jbvvvpvvvvvu\nlisS1uRi9+7dPPvsswB07doVf39/jhw5QlhYmENjdbZaHTfrbMSkDpWUlKguXbqojIwM9fPPP1sc\nuN6zZ88tO1hrTS5OnDihunbtqvbs2eOkKB3DmlxUFBsbq/797387MELHsSYXhw8fVkOHDlXXrl1T\nV65cUYGBgerQoUNOith+rMnF448/rhYuXKiUUio3N1f5+vqq8+fPOyNcu8vIyLBq4Nra46ZLtiSq\nu15ixYoVADz66KOMGTOGTz/9lG7dutGsWTPeeecdJ0dtH9bk4qWXXqKgoMDcD+/h4UFycrIzw7YL\na3LRUFiTi169ejFq1Cj69u2Lm5sbM2bMICAgwMmR1z1rcvHMM88wZcoUgoODKSsr4y9/+QutW7d2\ncuR1Lzo6mi+//JK8vDw6derEiy++SElJCVD742a9u+mQEEIIx3HJs5uEEEK4BikSQgghqiVFQggh\nRLWkSAghhKiWFIl6Jj4+noCAACZPnmzzthISEsjJyTGvz5gxg8OHD9u83eoUFxczbNgwQkJC+OCD\nD+y2H0cbNGiQxdfUZW6t2V9F1/+djUYj+fn5dRJLbX333Xd89tlnFl+XkpLCH/7wBwdEJKpVh6fn\nCgfo1auXys7OvuH52kxWFhERofbt21cXYVllz549atiwYQ7bX0WlpaVO2a8ruP7vbDQaVV5enhMj\nUuqdd95Rs2bNcug+G/J3wBZSJOqRRx99VDVu3FgFBQWppUuXqoULF6qHHnpIDRo0SD344IMqMzNT\nDR48WIWGhqrQ0FC1e/du83vj4uJUUFCQCg4OVvPmzVMbNmxQzZs3Vz179lQhISGqqKhIDRkyxHww\nee+991RQUJAKDAxUTz/9tHk7zZo1U88++6wKDg5W4eHh6syZMzfEef78efWb3/xG9e3bV4WHh6uD\nBw+qs2fPqq5duypPT09lMplUenp6pfckJyeroKAgZTKZ1B//+EfzxUAZGRlVfqbTp0+rwYMHK5PJ\npAIDA9XXX399Qxx33HGHevrpp1VoaKhKTExU//3vf9Wdd96pQkND1f33368uX75sft38+fOVyWRS\n/fr1UykpKWr48OGqa9eu6h//+IdSSqlLly6poUOHqtDQUBUUFKQ2bdpUKSdKKbVjxw41ZMgQ9dvf\n/lb16tVLTZo0yfyaIUOGqJSUlJvm8Mcff1QDBgxQQUFB6tlnn1XNmzev8ntgzf7KffDBBzf8nY1G\no3rhhRfMn+WHH35QSil1+fJlNWXKFNW/f38VEhJS6TNWdP13SSmlDhw4oAYMGKD69u2r7rvvPlVQ\nUGD+3E8//bTq37+/6tGjh/r6669VcXGx6tSpk2rXrp0ymUxq3bp1Kjk5Wd15550qJCREDRw4UB05\ncsT8GX/9618rpbSZXKdMmaIiIiJUly5dVHx8vDmmNWvWqP79+yuTyaQeffRRc0Fo1qyZevLJJ1Vw\ncLDatWtXlZ9H3JwUiXrGaDSarxR94YUXVFhYmPrpp5+UUkpdvXrV/Pjo0aMqLCxMKaXUp59+qgYO\nHKiKioqUUsr8DzgiIsJ84Kq4np2drTp37qzy8vLUtWvX1D333KM+/PBDpZRSBoNBffzxx0oppZ56\n6in18ssv3xDjrFmz1EsvvaSUUuqLL75QJpNJKaVUUlKS+R/89fr06aP27t2rlFJq3rx5Kigo6Kaf\n6bXXXlN//vOflVJKlZWVqUuXLlWZq1dffVUppU0Pfffdd6urV68qpbQDXXmMRqPRXAwef/xxFRQU\npC5fvqzOnTunfHx8lFLadNQXL140b6tbt27m/ZQfzHfs2KE8PT1Vdna2KisrU3feeaf5wFQx19Xl\n8N5771WJiYlKKaX+8Y9/VFskbra/nTt33vD66//ORqNRvfnmm0oppZYvX66mT5+ulFJq/vz56t13\n31VKad+RHj16qCtXrlTaVnXfpaCgIPXVV18ppZR6/vnn1dy5c837/uMf/2h+b3lLctWqVWr27Nnm\n7V68eFFdu3ZNKaXUtm3b1IQJE8yfsWKRGDRokCouLlZ5eXmqTZs26tq1ayotLU2NHTvW/P6ZM2eq\n1atXm3P9wQcfVJlHYR2XvOJaWMdgMBAVFUWTJk0Arc9/1qxZfPfdd7i7u3Ps2DEAPv/8c6ZOnUrT\npk0B8PLyMm9DXXctpVKKb7/9loiICNq0aQPApEmT+Oqrr/jNb35D48aNuffeewHo168f27ZtuyGu\nXbt2mW/2ExkZyfnz57l8+XK1k+0VFhZy+fJlBgwYAMCDDz7Ixx9/fNPP1L9/f6ZOnUpJSQnjxo0j\nODi4ym0/8MADAOzdu5e0tDQGDhxo3m75Y8A8t1NQUBBXrlyhWbNmNGvWjCZNmnDx4kVuu+025s+f\nz9dff42bmxunT5/m7NmzeF93p/n+/fvTsWNHAEwmE5mZmZX2A1Sbw71797J582ZAu3L2j3/8Y5Wf\nydL+qhqzuD7348ePByA0NNT8t9q6dSsfffQRr732GgA///wzWVlZ9OzZ0/y+7du33/BdunDhAhcu\nXGDw4MEAxMTEcP/991e5r8zMTHM8FWMqLCzk4Ycf5scff8RgMJivEq7IYDBw77334uHhQZs2bfD2\n9iY3N5ft27eTkpJinoepqKiI9u3bA+Du7s6ECRMs5lFUT4pEPXf77bebHy9dupQOHTqwZs0aSktL\nzf+Qy2dErUpVM0Be/5xSyvych4eH+Xk3NzeuXbtW5Xar2581Kr63us80ePBgvv76az7++GNiY2N5\n4oknqhzMb9asmfnx8OHDee+996rcZ3mhdXNzo3Hjxubn3dzcKCkp4T//+Q95eXns378fd3d3/P39\n+emnn6rdDmgHqKryY20OrWHN/uDGv2n5+65/z3/+8x+LEwBa+tte//vq9lXRggULGDp0KBs3buTE\niRPV3kWv4t+m4vZiYmJYtGjRDa9v2rTpLTk7tCPJ2U23kIsXL5r/B7V69WpKS0sB7eD4zjvvUFRU\nBEBBQQEALVq0MN+IpZzBYKB///58+eWXnD9/ntLSUhITExkyZIjVcQwePJi1a9cC2u0z27VrR/Pm\nzat9vZeXFy1atDDPN5WYmGj+h13dZzp58iTt2rVj+vTpTJ8+3eJtGAcMGMCuXbtIT08H4MqVK+ZW\nSUXVHQAvXryIt7c37u7u7NixgxMnTtx0f7URHh7Ohg0bAC0HdaWqv3NVRo4cSXx8vHm9qpxW9V3y\n9PSkVatW7Ny5E4A1a9ZYvFVqy5YtuXTpknn94sWL5hZRdfMJVfW3MRgMDB06lA0bNnDu3DkA8vPz\nOXny5E33L6wnRaKeuf5/RRXXH3vsMRISEjCZTBw5csR8YB45ciRRUVGEhYUREhLCkiVLAIiNjeX3\nv/89oaGhlf5X3L59e+Li4oiMjMRkMhEWFsbYsWNv2J/BYKjyf2kLFy4kJSWF4OBgnnnmGRISEm76\netBuljRjxgxCQkK4evUqLVu2vOln2rFjByaTidDQUNavX1/laZIV99WuXTtWrVpFdHQ0wcHBDBw4\nkCNHjlT5nqo+46RJk9i3bx99+/ZlzZo19O7du8r9WPO/1upy+Ne//pXXX38dk8lEenq6+W6Dlt5f\n3e/KVfd3vn7/CxYsoKSkhL59+xIYGMgLL7xww7aq+y4lJCTwpz/9ieDgYA4ePMjzzz9/09gjIyNJ\nS0sjJCSE9evX89RTTzF//nxCQ0MpLS2t8jNW9/3p3bs3L7/8MiNGjCA4OJgRI0aQm5tbbT5EzcgE\nf8IllI8DAMTFxXHmzBmWLl3q5Kgcq6ioiNtuuw3QWhLr1q1j48aNTo5KNHQyJiFcwieffMLixYu5\ndu0aRqORVatWOTskh0tJSWHWrFkopWjVqhVvv/22s0MSQloSQgghqidjEkIIIaolRUIIIUS1pEgI\nIYSolhQJIYQQ1ZIiIYQQolpSJIQQQlTr/wPtfJXPCprbqQAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x9cfc1cc>"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "If we increase the nozzle width as suggested by the blog post by a factor of 5, we manage to empty one canister in:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "one_canister_empty_time = time_to_empty(0.1, P_0, rho_0, V, A * 25) \n",
      "print '%.2f s' % one_canister_empty_time"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "4.88 s\n"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This is indeed an improvement, as the total emptying time now becomes"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print '%.2f s' % (number_canisters * one_canister_empty_time)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "53.68 s\n"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In the blog post however, the author suggests to use an apparatus to release all canisters at once. This yields a total emptying time equal to the one it takes to empty only one canister:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print '%.2f s' % one_canister_empty_time"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "4.88 s\n"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "However, by emptying the gas so quickly into the balloon, a drag force will appear. This drag force is usually expressed like:\n",
      "\n",
      "$$\n",
      "F_d = C_d \\; \\left( 1/2 \\; \\rho \\; u^2 \\; A \\right)\n",
      "$$\n",
      "\n",
      "where $F_d$ is the drag force, $C_d$ the drag coefficient, $\\rho$ the density of the fluid where the drag force is created, $u$ the speed at which the object is moving and $A$ the cross section in the flow direction of the object.\n",
      "\n",
      "Wikipedia has [a nice page about the drag coefficient](http://en.wikipedia.org/wiki/Drag_coefficient). See also [this page](http://wright.nasa.gov/airplane/drageq.html) from NASA. \n",
      "From Wikipedia we learn that a streamlined body has a drag coefficient of $0.04$ and that for a sphere it's $0.47$. Therefore, taking $0.1$ seems reasonable for a start."
     ]
    }
   ],
   "metadata": {}
  }
 ]
}