Problem Set 3 MECH 430

gistfile1.txt

{
 "metadata": {
  "name": "",
  "signature": "sha256:c7b95a8879669ca52a404e446df969830dd8b4ff13fdca2826dbffb744c886ae"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "IPython Notebook\n",
      "================\n",
      "\n",
      "If you're viewing this document as PDF, the input code will be hidden. If you desire to view it, you can do so at :bit.ly/selimb_PS3_MECH430\n",
      "\n",
      "Question 1\n",
      "=========="
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import sympy as sp\n",
      "sp.init_printing()\n",
      "from IPython.display import Math\n",
      "import matplotlib"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "A,P0,T0,gamma,R,M = sp.symbols('A P_0 T_0 gamma R M',real=True,positive=True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mdot = A*P0/sp.sqrt(T0)*sp.sqrt(gamma/R)*M/(1+(gamma-1)/2*M**2)**((gamma + 1)/(2*(gamma-1)))\n",
      "Math(\"\\dot{m}(M) = \" + sp.latex(mdot))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "$$\\dot{m}(M) = \\frac{A M P_{0} \\sqrt{\\gamma}}{\\sqrt{R} \\sqrt{T_{0}}} \\left(M^{2} \\left(\\frac{\\gamma}{2} - \\frac{1}{2}\\right) + 1\\right)^{- \\frac{\\gamma + 1}{2 \\gamma - 2}}$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 4,
       "text": [
        "<IPython.core.display.Math at 0x7bc78b0>"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "-----------\n",
      "\n",
      "To obtain the mach number corresponding to the *maximum* flow rate, we must solve the following equation for M :\n",
      "\\begin{equation}\n",
      "\\frac{\\text{d}\\left(\\dot{m}(M)\\right)}{\\text{d}M} = 0\n",
      "\\end{equation}\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dmdot_dm = sp.diff(mdot,M)\n",
      "Math(\"\\\\frac{\\\\text{d}\\\\left(\\\\dot{m}(M)\\\\right)}{\\\\text{d}M} = \"+ sp.latex(dmdot_dm))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "$$\\frac{\\text{d}\\left(\\dot{m}(M)\\right)}{\\text{d}M} = - \\frac{2 A M^{2} P_{0} \\sqrt{\\gamma} \\left(\\frac{\\gamma}{2} - \\frac{1}{2}\\right) \\left(\\gamma + 1\\right) \\left(M^{2} \\left(\\frac{\\gamma}{2} - \\frac{1}{2}\\right) + 1\\right)^{- \\frac{\\gamma + 1}{2 \\gamma - 2}}}{\\sqrt{R} \\sqrt{T_{0}} \\left(2 \\gamma - 2\\right) \\left(M^{2} \\left(\\frac{\\gamma}{2} - \\frac{1}{2}\\right) + 1\\right)} + \\frac{A P_{0} \\sqrt{\\gamma}}{\\sqrt{R} \\sqrt{T_{0}}} \\left(M^{2} \\left(\\frac{\\gamma}{2} - \\frac{1}{2}\\right) + 1\\right)^{- \\frac{\\gamma + 1}{2 \\gamma - 2}}$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 5,
       "text": [
        "<IPython.core.display.Math at 0x6016bf0>"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "----------------\n",
      "\n",
      "Solving for $M$:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Mmax = sp.solve(dmdot_dm,M)[0]\n",
      "print \"a)\"\n",
      "Math(\"M = \" + sp.latex(Mmax))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "a)\n"
       ]
      },
      {
       "latex": [
        "$$M = 1$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 6,
       "text": [
        "<IPython.core.display.Math at 0x7cfad10>"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "------------"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print \"b)\"\n",
      "Math(\"\\dot{m}(1) = \" + sp.latex(sp.simplify(mdot.subs(M,1))))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "b)\n"
       ]
      },
      {
       "latex": [
        "$$\\dot{m}(1) = \\frac{A P_{0} \\sqrt{\\gamma}}{\\sqrt{R} \\sqrt{T_{0}}} \\left(\\frac{\\gamma}{2} + \\frac{1}{2}\\right)^{- \\frac{\\gamma + 1}{2 \\gamma - 2}}$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 7,
       "text": [
        "<IPython.core.display.Math at 0x7c858d0>"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "-----\n",
      "\n",
      "Question 2\n",
      "=========="
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import sys\n",
      "sys.path.append('tools')\n",
      "import difftotal"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Assuming a *calorically perfect gas*, we derived:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "T, t, dT, dM, P, dP = sp.symbols('T t dT dM P dP')\n",
      "T0_eq = T*(1 + (gamma-1)/2*M**2)\n",
      "Math(\"T_0 = \" + sp.latex(T0_eq))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "$$T_0 = T \\left(M^{2} \\left(\\frac{\\gamma}{2} - \\frac{1}{2}\\right) + 1\\right)$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 9,
       "text": [
        "<IPython.core.display.Math at 0x7c77b30>"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "If the process is assumed to be adiabatic, $\\text{d}T_0 = 0$. Taking the total derivative of the above equation yields:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "diffT0 = difftotal.difftotal(T0_eq, t,\n",
      "                       {T:dT, M:dM})\n",
      "Math(\"\\\\text{d}T_0 = 0 = \" + sp.latex(sp.simplify(diffT0)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "$$\\text{d}T_0 = 0 = M T dM \\left(\\gamma - 1\\right) + \\frac{dT}{2} \\left(M^{2} \\left(\\gamma - 1\\right) + 2\\right)$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 10,
       "text": [
        "<IPython.core.display.Math at 0x205f330>"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Simplify and solve for $\\frac{\\text{d}T}{T}$:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dt_tM = sp.solve(diffT0,dT)[0]/T\n",
      "\n",
      "Math(r\"\\frac{\\text{d}T}{T} = \" + sp.latex(sp.simplify(dt_tM)) + \"~~~~(1)\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "$$\\frac{\\text{d}T}{T} = - \\frac{2 M dM \\left(\\gamma - 1\\right)}{M^{2} \\gamma - M^{2} + 2}~~~~(1)$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 11,
       "text": [
        "<IPython.core.display.Math at 0x7c6f4f0>"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We also derived in class the relationship between pressure change and temperature change for an isentropic process."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Math(r\"\\frac{\\text{d}T}{T} = \\frac{\\gamma - 1}{\\gamma}\\frac{\\text{d}P}{P} ~~~~(2)\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "$$\\frac{\\text{d}T}{T} = \\frac{\\gamma - 1}{\\gamma}\\frac{\\text{d}P}{P} ~~~~(2)$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 12,
       "text": [
        "<IPython.core.display.Math at 0x7c920f0>"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Equate $(1)$ and $(2)$ and obtain the following relationship between mach number change and pressure change:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dt_tP = (gamma-1)/(gamma)*(dP/P)\n",
      "M_P_eq = sp.Eq(dt_tM,dt_tP)\n",
      "Math(sp.latex(M_P_eq)+\"~~~~~~(*)\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "$$\\frac{2 M dM \\left(- \\gamma + 1\\right)}{M^{2} \\gamma - M^{2} + 2} = \\frac{dP \\left(\\gamma - 1\\right)}{P \\gamma}~~~~~~(*)$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 13,
       "text": [
        "<IPython.core.display.Math at 0x5b95910>"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "From Eq. (4.34) in [1], we can write a differential equation for the pressure as a function of the area.\n",
      "\n",
      "We can then substitute $dP$ in $(*)$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "rho, V, dA = sp.symbols('rho V dA')\n",
      "dP_eq = (rho*V**2*dA)/(A*(1-M**2))\n",
      "M_A_eq = M_P_eq.subs(dP,dP_eq)\n",
      "Math(\"\\\\text{d}P = \" + sp.latex(dP_eq) + \"~~\\\\Rightarrow ~~\" + sp.latex(M_A_eq))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "$$\\text{d}P = \\frac{V^{2} dA \\rho}{A \\left(- M^{2} + 1\\right)}~~\\Rightarrow ~~\\frac{2 M dM \\left(- \\gamma + 1\\right)}{M^{2} \\gamma - M^{2} + 2} = \\frac{V^{2} dA \\rho \\left(\\gamma - 1\\right)}{A P \\gamma \\left(- M^{2} + 1\\right)}$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 14,
       "text": [
        "<IPython.core.display.Math at 0x6016ef0>"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dA_A_temp = sp.solve(M_A_eq,dA)[0]/A\n",
      "Math(\"\\\\Rightarrow ~ \\\\frac{\\\\text{d}A}{A} = \" + sp.latex(sp.simplify(dA_A_temp)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "$$\\Rightarrow ~ \\frac{\\text{d}A}{A} = \\frac{2 M P dM \\gamma \\left(M^{2} - 1\\right)}{V^{2} \\rho \\left(M^{2} \\gamma - M^{2} + 2\\right)}$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 15,
       "text": [
        "<IPython.core.display.Math at 0x7ce8170>"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Finally, from ideal gas law we can easily prove that $\\rho V^2 = \\gamma \\cdot M^2 \\cdot P$. \n",
      "\n",
      "This allows us to simplify further and obtain a differential equation relating change in area to change in mach number."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dA_A = dA_A_temp.subs(rho*V**2,gamma*M**2*P);\n",
      "# dA_A\n",
      "#We don't print dA_A because sympy doesn't format the equation in a \n",
      "#nice and concise way. They are equivalent, however.\n",
      "Math(r\"\\frac{\\text{d}A}{A} = f(M)\\frac{\\text{d}M}{M}\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "$$\\frac{\\text{d}A}{A} = f(M)\\frac{\\text{d}M}{M}$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 16,
       "text": [
        "<IPython.core.display.Math at 0x7cfa750>"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Math(r\"~~~\\text{where},~~ f(M) = \\dfrac{M^2 - 1}{(1+\\frac{\\gamma-1}{2}M^2)}\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "$$~~~\\text{where},~~ f(M) = \\dfrac{M^2 - 1}{(1+\\frac{\\gamma-1}{2}M^2)}$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 17,
       "text": [
        "<IPython.core.display.Math at 0x6010a30>"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# sp.simplify(sp.integrate(dA_A/dM,M))\n",
      "Astar = sp.symbols('Astar')\n",
      "LHS = sp.exp(sp.integrate(1/A,(A,Astar,A)))\n",
      "RHS = sp.exp(sp.integrate(dA_A/dM,(M,1,M)))\n",
      "Math(r\"\\frac{A}{A_{star}} = \" + sp.latex(sp.simplify(sp.expand(RHS))) + \" = X(M)\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "$$\\frac{A}{A_{star}} = \\frac{1}{M} e^{\\frac{1}{2 \\gamma - 2} \\left(- \\gamma \\log{\\left (\\frac{\\gamma + 1}{\\gamma - 1} \\right )} + \\gamma \\log{\\left (\\frac{M^{2} \\left(\\gamma - 1\\right) + 2}{\\gamma - 1} \\right )} - \\log{\\left (\\frac{\\gamma + 1}{\\gamma - 1} \\right )} + \\log{\\left (\\frac{M^{2} \\left(\\gamma - 1\\right) + 2}{\\gamma - 1} \\right )}\\right)} = X(M)$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 18,
       "text": [
        "<IPython.core.display.Math at 0x7c8f0b0>"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "While this does not look like the required equation, we will prove that they are equal by checking values for two particular $M$ and $\\gamma = 1.4$. We will also plot both functions. \n",
      "\n",
      "The relation given in the question was:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "a = 2/(gamma+1); b = 1 + (gamma-1)/2*M**2; powr = (gamma+1)/(2*(gamma-1));\n",
      "RHS_desired = 1/M*(a*b)**powr\n",
      "Math(r\"\\frac{A}{A_{star}} = \" + sp.latex(RHS_desired) + \" = Y(M)\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "$$\\frac{A}{A_{star}} = \\frac{1}{M} \\left(\\frac{1}{\\gamma + 1} \\left(2 M^{2} \\left(\\frac{\\gamma}{2} - \\frac{1}{2}\\right) + 2\\right)\\right)^{\\frac{\\gamma + 1}{2 \\gamma - 2}} = Y(M)$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 19,
       "text": [
        "<IPython.core.display.Math at 0x7cfadf0>"
       ]
      }
     ],
     "prompt_number": 19
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Math(r\"X(2) = \" + sp.latex(RHS.subs({gamma:1.4,M:2}).evalf()) +  \",~ X(1) = \" + sp.latex(RHS.subs({gamma:1.4,M:1}).evalf()))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "$$X(2) = 1.6875,~ X(1) = 1.0$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 20,
       "text": [
        "<IPython.core.display.Math at 0x7e88ab0>"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 21
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Math(r\"Y(2) = \" + sp.latex(RHS_desired.subs({gamma:1.4,M:2}))  +  \",~ Y(1) = \" + sp.latex(RHS_desired.subs({gamma:1.4,M:1}).evalf()))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "$$Y(2) = 1.6875,~ Y(1) = 1.0$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 22,
       "text": [
        "<IPython.core.display.Math at 0x853fd90>"
       ]
      }
     ],
     "prompt_number": 22
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "rhs_desired_plot = sp.plotting.plot(RHS_desired.subs(gamma,1.4),\n",
      "                                    (M,0.5,2),\n",
      "                                    axis_center = (0.5,1),\n",
      "                                    title='Y(m)',\n",
      "                                    ylabel='A/A_star',\n",
      "                                    line_color = 1000)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEfCAYAAABF6WFuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4VGXexvFv6L1EikgxUgQCUkSaoAxCKCIgRaRIpIiA\nUtRV1HV3yWtDWKQJAoIBAQVUJIDSm0aqIFUIRUCKEJp0CCnz/nGCZmPKJJmZZ8r9ua65mHIyc4s5\n/OapB0RERERERERERERERERERERERERERERERERERDxUduBl0yFERMSz/GQ6gIiIuEd2B48LBtoC\nV4FiQKnE22kX5RIREUMCHDxuPWBP4fmmzosiIiIi4l2CyXy3aklgH5DLeXFEXMPRFgPAE1gnRp4k\nz73t3DgiRs0BbgN9kjzXBFgAVAcmAfOBLzP5/pOA/cDELGQU8RhTgVnASWA4sBf41GgiEecLxBo3\na574OA9wEAjFGlO7QNa+8T8M7MlKQBFPcueXeXfinwWAHw1lEXGlzsARIB8wAvgu8flQYGWyY48B\nr2KdF1exviyVBJYBl4FVQJEkx+cArgNlXRNdxL22Jv65GSiN9U3qsLk4Ii71NbAYOI/1+w7wX+Cj\nZMcdBTYCxYF7gGjgZ6AmkBtYA/wn2c/swprhJ+Kxcjh43LdAUayTY3vic9NckkjEvBeAX4F/AqcS\nnyuM1ZWU3EfAucT7kVjFYVfi44VAs2THX018LxGP5WhhGAXcwhqE+w6rxXDLVaFEDDuL1Vr4Jclz\nfwAFUzg2Osn9m8ke38Lqdk2qIHDJCRlFXCabg8dtTHL/FtYv9sZUjhXxRbuB+x04Lq2ZfjmAivzV\nohDxSOm1GEph9Z3mAx7E+qW3A4USnxPxF6uB8Vizkm5n8j3qYQ1Yn3BSJhGXSK8wtAB6YQ3AfZjk\n+atY/a8i/iIaWAs8SdrrGOzJ7id93AOY7PxoImZ0yuTPhWOdUKnN3X4V2JF42wPE8b/T+0Q8SVX+\nmqGXUSXQymfxMS9hdR8FYM3V/hlo6cDPPQLUxrFFPU9gNddFRMQL3FnY1hJrCl51rG/5jgjCscLw\nBdA3w8lERMSpHJ2VdGemRRtgNtaWGM6UD6voLHDy+4qISAY5Whi2Y20H8DiwHKtbKcGJOdpibbGh\n+d0iIoY5usCtL1ALazXoDeAuoHeS16vxv4uBMqorMDe1F5999ll7UFDQn49tNhs2my0LHyci4vsC\nAgIysoP2Xz/npM/fgTXInJIgYAnwQCqvF8batKwM1srRlNjt9pSuEyQiIqnJbGFwtMWQWXOx9rMv\nhrWoZziQM/G1qYl/PgmsIPWiICIibuSOFoMzqMUgIpIBq1b9SosWFTP1b7yjg88iIuIl4uISeOGF\npZn+eWcVhhgnvY+IiGTRnDm7KV06pc2AHZOZZkZFoBvWTKJqmf7kjFFXkoiIA2Jj46lSZRLh4e2w\n2e5zaVdSaeAV4CesxW3ZsQqDiIh4kFmzdhEUVIQmTYIy/R7pVZP+WK2DEliXO/wK65KH92X6EzNH\nLQYRkXTcvh1P5coTmT27A40bl3PZdNWJWCudh6KLi4iIeLSZM3dSqVIgjRuXy9L7OHKhnqeACfzV\nasiZ5k+IiIjbxcTE8d57kcybl9mrJPwlvTGG81gXFmmCddGey1jXV4gC3s/yp4uIiFOEh+8gOLg4\nDRuWzfJ7pdf/dA/wewrP3481+Px2lhM4RmMMIiKpiImJo2LFj1iwoAv16pX+83lXjTF8CgQC67DG\nGn7EusraQdxXFEREJA3Tp/9MjRol/6coZEV6haE1kBewAR2B0Vh7Hi3DKhTHnZJCREQy5ebNWFas\n+JWwsCZOe09HNtG7iVUIliU+Lg88A0zB2hyvntPSiIhIhkyevI3s2bNRt65zWguQsd1VH8Ra09AF\nOIp1tbWpaf6EiIi4zNWrMYwcuYE1a0Kd+r7pFYbKWMXgaeAc1gK3AKyuJRERMWjcuM00b16e6tVL\nOPV90xuxTgC+BQbx13jCUbTyWUTEqIsXb3L//R+xefNzVKwYmOIxmZ2VlN46ho5YYww/YI0pNMN5\n13AQEZFM+u9/N9CxY9VUi0JWOPqPfAGgPVa3UlNgFrAQWOn0RClTi0FEJNGZM9eoVu1jdu7sT9my\nhVM9zp3XfA4EOmMtcHssMx+aCSoMIiKJhg5dRkBAAOPGtUrzOHcWBhNUGEREgOPHL1O79lT27XuB\nkiULpHmsq8YYRETEg0yatJX+/eukWxSyIiPrGERExKCoqPOEh+/kwIEXXfo5ajGIiHiJt95ay6uv\nNiQwMJ9LP0ctBhERL7B16ym2bDnJ7NkdXP5Zrm4xhGNdv2FPGsfYgB1Y15Je7+I8IiJex26388Yb\nqxk+vAn58rn+WmmuLgwzgLTmUxUBJgFtgepY02BFRCSJlSt/5dSpq/TuXdstn+fqwhAJ/JHG692x\nNuM7mfj4vIvziIh4lYQEO2+8sYb33nuMHDncMyxsevC5En9dCGgb0NNsHBERzzJ//l5y5sxGp05V\n3faZpgefc2Jt590MyAdsAjYDh5IfGBYW9ud9m82GzWZzS0AREVNu347nX/9ax7RpbcnkWrVMMV0Y\nTmB1H93kr836apJOYRAR8Qdz5uymYsVAHnvMvRtam+5KWgQ0BrJjtRjqA/uMJhIR8QBXr8bw1ltr\nGTmymds/29UthrlAE6xLgJ4AhmN1H4F19bcorGtH78a69sM0VBhERBg1agMtWlSgVq1Sbv9sbaIn\nIuJhTp26Qo0aU9ixoz/lyqW+rXZ6tLuqiIiPeO65xRQrlo8PPmiepffJbGEwPfgsIiJJ7N17liVL\nDnLw4CBjGUwPPouISBLDhq3in/9sTOHCeYxlUItBRMRDrFlzhIMHLxAR0dVoDrUYREQ8QEKCndde\nW8WIEc3IlSu70SwqDCIiHuCLL/aQK1d2OncONh1FXUkiIqbduBHLiBGRTJ/ezq1bX6RGLQYREcPG\njNlEcHAJGjYsazoKoBaDiIhRv/9+lbFjN/PTT/1MR/mT+TaLY7TATUR8Up8+iyhePB8jR4Y4/b21\nwE1ExMv8/PNpli07zIED5hazpURjDCIiBtjtdl55ZQVhYU0oVCi36Tj/Q4VBRMSAiIgoLly4Sd++\nD5qO8jfqShIRcbOYmDhee20Vkye3cdt1nDPC8xKJiPi4SZN+okqVYoSEVDAdJUVqMYiIuNHZs9dZ\nteoIY8a0MB0lVWoxiIi40VtvraFKlbuoWrW46SipUotBRMRNtm//nSVLDhIV5VnTU5NTi0FExA3s\ndjtDhy7nnXeaUqSIuWstOEKFQUTEDebO3cuNG7H06VPbdJR0qStJRMTFrl+/zeuvr2bu3E5kz+75\n38c9P6GIiJcbMeJHHn30Xho3Lmc6ikPUYhARcaEjR/5gypRt7Nw5wHQUh7m6xRAORAN7UnndBlwG\ndiTe/uXiPCIibhUWtp7//KcJZcoUMh3FYa4uDDOAVukc8z1QO/H2rovziIi4zbJlh9i06ST9+9cx\nHSVDXF0YIoE/0jnGW64JISLisJiYOIYMWc6ECa3Indu7eu1NDz7bgYeBXcBSwPxVsEVEnGD06I1U\nq1ac1q0rmY6SYabL2M9AWeAG0BqIAO5P6cCwsLA/79tsNmw2m+vTiYhkwm+/XWLs2M1s2/a86SiZ\n4o5unCBgCfCAA8ceBeoAF5M9r0t7iojX6NTpS2rWLMl//tPEaA5vvbRnSeAsVpdSPaxClbwoiIh4\njTVrjrBv3znmzOlgOkqmubowzAWaAMWAE8BwIGfia1OBzsBAIA6rO6mri/OIiLjMrVtxDBjwHePG\ntSRv3pzp/4CH8pYZQepKEhGP9/bb37NjxxkWLnzadBTAe7uSRER8wuHDF5kwYQs//9zfdJQsMz1d\nVUTE69ntdgYNWsrrrzeiXLnCpuNkmQqDiEgWff31Pk6evMJLLzUwHcUp1JUkIpIFV67E8PLLK5g7\ntxM5c2Y3Hccp1GIQEcmCsWM3ERJSgUceudd0FKdRi0FEJJO2b/+dSZN+Yu/egaajOJVaDCIimRAX\nl8Dzz3/LqFEhlChRwHQcp1JhEBHJhI8+2kLhwrl59tmapqM4nRa4iYhk0PHjl3nwwals3NiX+++/\ny3ScVGV2gZtaDCIiGWC323nxxaUMHVrfo4tCVmjwWUQkAxYs2M+vv15kwYIupqO4jAqDiIiDLl26\nydChy5k/vzO5cvnGmoWUaIxBRMRBAwd+S7Fi+XjnncdMR3GINtETEXGh778/xpIlB9m79wXTUVxO\ng88iIum4eTOW555bwqRJj1OkSB7TcVxOhUFEJB3/93/f8+CDpWjfvorpKG6hriQRkTT8/PNpZszY\nye7dA0xHcRu1GEREUhEbG0/fvov5739DKFnSt7a9SIsKg4hIKkaP3kjJkvnp2bOG6ShupcIgIpKC\n/fvPERl5nKlTnyCTsz69lgqDiEgy8fEJ9O69iCeeuJ977y1iOo7bqTCIiCQzZswm8ubNyYABD5mO\nYoRmJYmIJBEVdZ6RIzfw00/9yJbNv7qQ7nB1iyEciAb2pHNcXSAO6OjiPCIiqbrThfR//2fjvvuK\nmo5jjKsLwwygVTrHZAdGAsvxnr2bRMQHjRu3mTx5cjBwYF3TUYxydWGIBP5I55jBwNfAORdnERFJ\n1YED55k8eRufftrOb7uQ7jA9+FwaaA9MTnysLVRFxO3i4hIIDY3gH/9oSPny/tuFdIfpwedxwBtY\nBSGANLqSwsLC/rxvs9mw2WwujiYi/mLUqA0UKpSb/v39cxZScu5oLwUBS4AHUnjtSJIMxYAbQD9g\ncbLjdD0GEXGJXbvO0Lz5bH7++XnKli1sOo5Teev1GMonuT8Dq4AkLwoiIi4RExNHaGgEo0eH+FxR\nyApXF4a5QBOs1sAJYDiQM/G1qS7+bBGRNH3wwY8EBRUhNLSm6SgexVuG3tWVJCJOtWHDcbp3/4Yt\nW57j7rt9c+fUzHYlmZ6VJCLidlevxhAaGsGECa18tihkhVoMIuJ3nnvOGsqcPr2d4SSu5a2DzyIi\nbhUREcW6dcfYubO/6SgeS4VBRPzGmTPXGDDgWxYs6ELBgrlNx/FY6koSEb9gt9tp0+YLbLZ7GTas\nsek4bqHBZxGRNEya9BPnz9/g5Zcbmo7i8dRiEBGft3fvWZo2/YyNG/tQqdJdpuO4jVoMIiIpuHUr\nju7dFzByZHO/KgpZocIgIj7tgw9+pHbtu+ndu5bpKF5Ds5JExGctWXKAGTN2smNHfzLZq+KXVBhE\nxCedOnWFfv2W8PXXXQgMzGs6jldRV5KI+Jz4+ASeeWYhL75Yl8aNy5mO43VUGETE54wY8SN2u51/\n/vMR01G8krqSRMSnbNhwnI8+2sr27c+TPbu++2aG1/ytxcUlmI4gIh7u4sWbvPdeJNOmtaVMmUKm\n43gtrykMy5cfNh1BRDyY3W6nV68IqlQpRrt2lU3H8WpeUxjCw3eYjiAiHmzMmE1ER1/ngw+am47i\n9bymMKxde5To6GumY4iIB9q06QSjRm1k/vzO5MqV3XQcr+c1heHJJ6swZ85u0zFExMNcvHiTrl0X\nMG1aW4KCipiO4xO8pjD06VOb8PCdaDM9EbkjIcHOW2+toXPnqhpXcCKvKQyPPFKO27fj2br1lOko\nIuIhRo3awK5d0YwYoXEFZ/KawhAQEEDv3rU0CC0igDXuOH78Fr788imNKziZqwtDOBAN7Enl9fbA\nLmAHsB14LK03e/bZmnz11T5u3Ih1akgR8S6nTl2hR49vmDOng9YruICrC8MMoFUar68GagK1gV7A\nJ2m9WenShWjYsCwLFuxzWkAR8S6xsfF06fI1gwfXo1mz8qbj+CRXF4ZI4I80Xr+e5H4B4Hx6b9in\nTy3Cw3dmNZeIeKlhw1YRGJiXN97wj+s2m+AJYwxPAvuBZcCQ9A5u27YyefPm4ODBdGuIiPiYiIgo\nli07zKxZT5Itm66v4CqeUBgigKpAW2B2egfnypWdGjVKMmXKdpcHExHPsWdPNP36LWHevE4ULarr\nK7iSJ+2uGomV5y7gQvIXw8LC/rxfq9ZDvPjiPt55pyn58+dyW0ARMePixZt06DCf8eNbUatWKdNx\nfJ472mJBwBLggRReqwAcAezAg8BXic8lZ0++sO3JJ+fx+OOVeP75Ok4NKyKeJT4+gSeemEvVqsUY\nM6al6TheJSCT1zN1dVfSXGAjUBk4AfQB+ifeADphTWXdAYwHujr6xoMG1WPixK1aCS3i4/7zn3Xc\nuhXHqFEhpqP4DW8Zvflbi8Fut1O16iSmTWvLI4/cayiWiLjSggX7eOWVlWzb1o/ixfObjuN1PLXF\n4DIBAQGJrYafTEcRERfYvTua996L5JtvuqgouJnXFgaA0NCarFr1K6dOXTEdRUSc6Ny567RvP4/X\nXnuYOnXuMR3H73h1YShUKDfdulXnk080dVXEV9y+HU/nzl/RrVt1unVLac6KuJrXjjHcsW/fOZo1\nm8Vvv72kjbREvJzdbmfAgG85c+Y6Cxc+rUVsWeR3Ywx3BAcXJzi4OF9/rf2TRLzdxx//xIYNJ5gz\np4OKgkFeXxgA/vGPBowdu1lTV0W82Lp1R3nnnR9YvLgbBQvmNh3Hr/lEYWjVqhLXr99m7dqjpqOI\nSCZERZ0nNDSCr756ivLli5qO4/d8ojBkyxbAsGGN+OCDDaajiEgGnT9/gyee+IK337ZpTZKH8InC\nANC9+wNERZ1n+/bfTUcREQfFxMTRocN8nnoqmN69a5uOI4m8ZXQn1VlJSY0du4lNm07y5ZdPuSGS\niGSF3W4nNDSCmzdj+fLLpzTY7AJ+OyspqX796rBu3TEOHfrb5qwi4mHeffcHDhw4z6xZmoHkaXyq\nMBQokIuBAx9i9OiNpqOISBo+/3w327efZvHibuTLl9N0HEnGW8q0Q11JYC2lr1x5Ir/88gKlShV0\ncSwRyai1a4/SrdsC1q4NpVq1Eqbj+DR1JSUqXjw/PXo8wPjxW0xHEZFk9uyJpmvXr5k/v7OKggfz\nuRYDwLFjl6hT5xOOHBlC4cJ5XBhLRBx18uQVHn74U0aNCqFr1+qm4/gFtRiSCAoqQseOVZgzZ7fp\nKCICXL58i9atP2fw4HoqCl7AJ1sMAPv3n6NJk5kcPjyEQoW0vF7ElJiYONq3n0elSoFMmNCaTH6J\nlUxQiyGZqlWL07JlRcaN22w6iojfSkiw8+yzEdx3XxHGjWulouAlvOX/UoZbDACHD1+kQYPpHDw4\nmMDAvC6IJSKpsdvtDBmyjD17zrJ8+TPkyZPDdCS/oxZDCipWDKRDhypa1yBiwLvv/kBk5HEWLeqq\nouBlfLrFAHD8+GVq157K/v0vUqKErhsr4g5Tpmxj9OiN/PhjH+6+u4DpOH5LLYZUlCtXmO7dqzNy\n5I+mo4j4hUWLDvDOOz+wYsUzKgpeyudbDACnT1+lWrWP2bNnIKVLF3JiLBFJaunSQ/TqFcGaNaE8\n8EBJ03H8nqe2GMKBaGBPKq/3AHYBu4ENQA1XhChVqiB9+9bm/fcjXfH2IgKsX3+MXr0iWLy4m4qC\nl3N1YZgBtErj9SPAo1gF4R3gE1cFGTasEStXHiEq6ryrPkLEb23ZcpIuXb5i3rzONGhQxnQcySJX\nF4ZI4I80Xt8EXE68vwVw2W9U8eL5GTCgDq++utJVHyHil3bvjqZdu3nMmNGexx67z3QccQJPGnzu\nCyx15QcMHlyfAwcusGLFYVd+jIjfOHDgPK1azeGjj1rTps39puOIk3hKYWgK9AFed+WH5MqVndGj\nQ3jllZXExSW48qNEfN7Ro3/Qp89i3nvvMbp0qWY6jjiRO2YlBQFLgAdSeb0G8A3WWERqX+Xtw4cP\n//OBzWbDZrNlKozdbqd589l06lSVF16om6n3EPF3x45dwmabyRtvNGLAAJ1Hniqzs5JMF4ZywFrg\nGSCtTY2yNF01ud27owkJmU1U1IsULaqtMkQy4tixSzRt+hmvvtqQF1+sZzqOpMFTC8NcoAlQDGva\n6nDgznX8pgLTgQ7A8cTnYoGUftOcWhgA+vdfQv78uRgzpqVT31fEl/32m1UUXn65AYMH1zcdR9Lh\nqYXBWZxeGM6evU61ah+zYUMf7r//Lqe+t4gvOn78MjbbTIYOrc/QoQ1MxxEHeOoCN49VokR+Xn+9\nkaavijjgxInLNG36GUOGqCj4A78tDACDB9cjIACWLDlgOoqIxzp27BJ9+y5m8OB6vPSSioI/8OvC\nkDt3DoYObcCgQcu4ejXGdBwRj3Po0AWaNJlJu3aVVRT8iN+OMSTVp88iChbMxfjxrV32GSLeZt++\nc7RoMZuwMBvPPfeg6TiSCRp8zoILF25QvfpkFi3qSr16pV32OSLeYufOM7Ru/TmjR4fQo4dL9rYU\nN9DgcxbcdVc+PvywBf36LSE2Nt50HBGjtm49RcuWc5g4sbWKgp9SYUjUrVt1SpUqwJgxm0xHETFm\n7dqjhIWtJzy8HZ06BZuOI4aoKymJo0f/oG7daWze/BwVKwa6/PNEPMmCBfsYOPA7vvzyKWy2INNx\nxAk0xuAkY8duYu/ec0yb1pZs2bzlr0ckaz75ZDthYev57rvu1K5dynQccRIVBieJj0+gSZOZdOxY\nlVdeaeiWzxQxxW638/77kYSH72TFimfUUvYxKgxOdPToH9SrN53Vq3tSs+bdbvtcEXdKSLDz8svL\nWb/+N5Yv70GpUgVNRxIn06wkJ7rvvqJ8+GELevT4hlu34kzHEXG627fj6NlzITt2nOH773upKMj/\nUGFIRc+eNQgOLs6bb642HUXEqS5fvsXjj39BkSJ5WLHiGYoUyWM6kngYFYZUBAQEMGXKE3z99X5W\nrvzVdBwRpzh+/DKNGoVTtWoxJkxoRd68OdP/IfE7KgxpCAzMy4wZ7enTZxEXLtwwHUckS3bsOM3D\nD39Knz61mTChNdmz6/SXlGnw2QGvvLKCGzduM3nyE2RyLEfEqKVLD9GrVwSTJ7fRwjU/ollJLnTr\nVhw2mzWFddiwRsZyiGSU3W5n4sStrFhxmLfeepSGDcuajiRupMLgYidOXKZevenMmvUkISEVjGYR\ncURsbDyDBy/jxx+Ps2RJN+67r6jpSOJmmq7qYmXLFmbu3E707LmQo0f/MB1HJE0XLtygRYs5nDp1\nlY0b+6ooSIaoMGSAzRbEG280pmPHL7lxI9Z0HJEU7d9/jvr1p1O37j1ERDxNoUK5TUcSL6OupAyy\n2+307LkQgNmzO2gwWjzKokVRTJ26jS5dqtOrVy3TccQwjTG40Y0bsTRqFE6vXjV1YXTxCPHxCQwf\nvp5Zs3bx1VdPUb9+GdORxANktjDkcHYQf5AvX06++aYLDRt+Ss2ad2uLYjHq4sWbdO++gJiYeLZt\ne54SJfKbjiReztVjDOFANLAnlderAJuAW8A/XJzFqe67ryizZ3fgww83cujQBdNxxE/t3HmGhx76\nhOrVS7BqVU8VBXEKVxeGGUCrNF6/AAwGRrs4h0uEhFSgXbvKtGr1OWfOXDMdR/zMvHl7CQmZzYgR\nzRg9ugU5cmguiTiHq7uSIoGgNF4/l3hr4+IcLtOvXx1On75G69af8/33vTQDRFzu+vXbDBq0jL17\nz7JuXSjVq5c0HUl8jL5iOMG///0oDRqUpkOH+cTEaJtucZ3du6N56KFp2O121q17VkVBXEKFwQkC\nAgKYOPFxihTJQ2hoBAkJnjODSnyD3W7nk0+206zZLN58szEzZz5JgQK5TMcSH+U1s5LCwsL+vG+z\n2bDZbMaypCR79mx8/nlHWracw0svLWf8+FZa4yBOcfnyLZ5//luios4TGdmbKlWKmY4kPs4d/3IF\nAUuAB9I4Jgy4CnyYyusetY4hLZcu3aJHj2+oW/cehg9vouIgWbJhw3FGjdpA6dKF+PDDFrp+gmSI\npy5wmws0AYphTVsdDtz5zZ4K3A38BBQCErCKQzCQfIqP1xQGgLNnrxMSMpuWLSswcmRzFQfJsNu3\n4xk+fB0zZ+5iypQ2tG9fxXQk8UKeWhicxasKA1ibmLVsOYdGjcoybpy6lcRxe/eepWfPhZQrV5hp\n09pqbYJkmnZX9TB33ZWP1atD2br1dwYO/E4D0pKuhAQ7Y8ZsomnTzxg0qC4REU+rKIgR3vI11uta\nDHdcvRpDmzZfUKFCINOnt9XlFCVFBw6c57nnllCxYiD//vejlC+vbbIl69Ri8FAFC+Zm2bIeHD9+\nmdDQCOLiEkxHEg8SF5fAyJE/0rjxDJ5+uhqfftpORUGMU4vBTW7ejKVfvyXcvh3P9OnttEJa2LXr\nDH37LqZo0bxMm9aWoKAipiOJj9HgsxeIjY1n0KClbNhwQpda9GPXrt3m7be/Z8+eaDp3DqZPn9qa\nnCAuoa4kL5AzZ3amTHmC55+vw8MPh/Pjj8dNRxI3stvtLFy4n+DgSZw5c42ZM5+kb98HVRTE43jL\nb6RPtBiSWr78MKGhCxk1KkRX2vIDR478weDByzh27BIff/w4TZoEmY4kfkBdSV5o375ztG07l86d\nq/L++800Y8kHXb9+m2nTfubdd39g2LBGvPRSA3Llym46lvgJdSV5oeDg4mzZ8hxbtpxiwIBvuXjx\npulI4iQJCXY++2wnVapMYvv239m+/XmGDWukoiBeQYXBsGLF8rFypXXlrZo1p7Bq1a+mI0kWrVt3\nlIce+oQpU7Yzf35nZs/uyL33asaReA91JXmQ1auP0Lv3Ijp1qsqIEc20YZqXOXDgPMOGrWb37mg+\n+KAZXbpU08CyGKUxBh9x8eJNBg78jr17z/L55x2pVetu05EkHcePX+a9935g06aT9OxZg8GD65Mn\nj9fsaC8+TGMMPiIwMC/z5nXizTcbExIym1GjNhAfr9XSnujEicsMHPgttWtPJTAwL2vXPstrrzVS\nURCvpxaDB/vtt0uEhkYQEADTprWlUqW7TEcS4OTJK4wYEcncuXvp1+9BXn31YYoX12Z34nnUYvBB\n995bhLVrQ+nWrToNGnzKv/61luvXb5uO5beOH7/EkCHLqFlzCvnz5+LAgUGMHBmioiA+R4XBw2XP\nno3+/R93N+eWAAAITElEQVRi164BHD16iSpVJjFv3l78sQVlyvbtv9O9+wJq1ZpKqVIF2L//RUaN\nUkEQ36WuJC8TGfkbQ4Ysp2DBXEyY0IpatUqZjuSTEhLsLF16iA8/3MThwxcZOrQ+/fo9SOHCeUxH\nE3GYZiX5kfj4BKZP/5nFiw9SpEge/v3vR3WBeCe5fPkWs2fvJiIiiosXb/Lqqw/z1FPB5MyphWni\nfVQY/NDVqzFMnLiVsWM3ExJSQQUiC7Zt+50pU7axYMF+WrSowODBdWnUqJzWIYhXU2HwY0kLRI8e\nD9C5czAPP1xW/6il48yZa8ybt5fNm0+yZcsp+vevQ+/etShZsoDpaCJOocIgXL0aw6xZuxg/fgv5\n8uXkhRfq0qPHA+TPn8t0NI9x40YsixZFMXv2bjZuPEG7dpUJDa3BY4+VJ1s2bzkdRByjwiB/Skiw\ns2bNESZN+onIyOP07FmDAQPqUKVKcdPRjLhyJYalSw+xcGEUly7dBALo2bMGHTpUUdEUn6bCICk6\nfvwyU6duY9euaKKjr9OlSzCdOwf7/NXjTp68wtq1R5k//xciI3+jceNydOxYlXbtKlOihKaZin/w\n1MIQDrQBzgIPpHLMBKA1cAPoBexI4RgVhiyKi0tg/fpjfPXVL3zzTRRBQUV46qlg2rW7n8qVi3n9\neMSVKzGsX3+M1auPsGrVEc6du06vXjWpU+ceHn+8kqaZil/y1MLwCHANmEXKheFxYFDin/WB8UCD\nFI5zWWFYv349NpvNJe/tCs7Ie6dILFiwjzVrjnLt2m1stqA/b5UqBTq1UDj779hut3PixBW2bj3F\n1q2nOHXqKosWRVG/fhlCQsoTElKe2rVLZWnMwB9/L9xNmV0vICCgKbA+oz/n6t2+IoGgNF5vB3yW\neH8LUAQoCUS7NtZfvO1/tDPy5siRjebNy9O8eXnsdjtHj15i/fpjrF9/jHfe+YH4+ASaNLmX+vXL\nEBxcnODg4pQuXTDTxSIrmWNj4zl69BIHD55n587oP4sBQP36ZahX7x7atKnEtGltyZfPeduU++Pv\nhbsps1vY8MDCkJ7SwIkkj08CZXBjYfB3AQEBlC9flPLli9KnT23sdjtHjvzBhg0n2LLlJIsXH2Df\nvnPcuBFL1apWkahatRhVqhSjePF8lCiRnxIl8lOgQK5MFY7Y2Hiio69z+vRVTp++RnT0NQ4dusiB\nAxc4cOA8x45d4p57ChIcXJxq1Yrz7LM1mTTpccqUKeT13V8insp0YYC/d2dpMMGggIAAKlQIpEKF\nQEJDa/75/MWLN9m//xz79p1j//7zHDx4gb17z3L27HXOnr1OfLydSpUCSUiwkzNndnLkyEbOnNnI\nkSMbJ07s4tSpJURFnefatdtcu3ab69etP4OCinD+/A1KlSpIqVIFuPfewpQuXYjQ0BpUrlyMihUD\ntY21iA8KAvak8toUoGuSx1FYXUnJ7cQqGLrppptuujl+m4mHCiL1wvA4sDTxfgNgszsCiYiIOXOB\n34HbWGMJfYD+ibc7JgKHgV3Ag+4OKCIiIiIiQius8YtDwOupHGPDWly3l0xM73KB9DIXA5Zjjb/s\nxVocaFI41myy1LoNwVrMeAirdVjbHaHSkV7mHlhZdwMbgBpuypUaR/6OAeoCcUBHlydKnyOZbXjO\nuZdeXk877wDKAuuAX7AyDUnlOE87/4zKjtVVFQTkxPofWjXZMUWw/lLLJD42vXe1I5nDgBGJ94sB\nFzA7y+wRrF82R8aT6uMZ40npZW4IFE683wrzmdPLC9bvzlrgW6CTO0KlI73MnnbupZc3DM867wDu\nBmol3i8AHODv/15k6Pzzh0t71sP6R/YYEAvMA9onO6Y7sABrHQXAeXeFS4UjmU8DhRLvF8L6BY1z\nU76URAJ/pPF6aosZTUov8ybgcuL9Lfz1j5cp6eUFGAx8DZxzfRyHpJfZ08699PJ62nkHcAbryyNY\nO03sB+5JdkyGzj9/KAwpLaIrneyYSkAgVnNsG9DTPdFS5UjmaUA1rMH9XcBQ90TLtNQWM3qLvvz1\njctTlcb6AjE58bHdYBZHedq5lx5PP++CsFo8W5I9n6Hzz3QTyB0cOTlyYs2Iagbkw/qmuBmrP84E\nRzL/E+tbgg2oAKwCagJXXRcry7x1MWNTrBl1jUwHScc44A2sv9cAvGP3ZE8799LjyeddAazW4lCs\nlkNyDp9//tBiOIU1OHNHWf5qtt5xAlgJ3MRqGv6A9T/bFEcyPwx8lXj/V+AoUNn10TIt+X9TmcTn\nPF0NrG+J7Ui/G8e0Oljdjkexxhc+xsrtyTzt3EuPp553ObG65OYAESm87q3nn8vkwPofGATkIuWB\n3CrAaqyBu3xYA0/B7ov4N45kHgMMT7xfEqtwBLopX2qC8L7FjEGknrkc1lhPSjv+mhJE+rOSAGbg\nGbOSIO3MnnbuQdp5PfG8C8DawXpsGsd46vlnVGuskfrDwJuJzyVfaPcq1uyIPaQ+3cud0stcDFiC\n1c+5B2sQzyRvXMyYXubpWN9idyTethrImJQjf8d3eEphcCSzJ5176eX1tPMOoDGQgPUF8s7vams8\n//wTERERERERERERERERERERERERERERERHvlwDMTvI4B9amikscfQN/2BJDRMSfXMfa6C9P4uMQ\nrBXaDu9NpsIgIuJ7lgJtEu93w1rR7fCmiioMIiK+Zz7QFcgNPMDft+FOkwqDiIjv2YO1GWA34LuM\n/rA/XI9BRMQfLQZGA02A4hn5QRUGERHfFI51DZFfsC4s5DB1JYmI+JY7s49OYW21fec5b7liooiI\niIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiI+Kv/B6wCZKWug0MuAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x8543650>"
       ]
      }
     ],
     "prompt_number": 23
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "rhs_plot = sp.plotting.plot(RHS.subs(gamma,1.4),(M,0.5,2),\n",
      "                            axis_center = (0.5,1),\n",
      "                            title='X(m)',\n",
      "                            ylabel='A/A_star')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEfCAYAAABF6WFuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X2cjPX+x/EXIklu4hQhqyI36U5HkjK6FaU7OrojJd2S\nOpXqdFpHSYUoR6XjJt1R8qOco1RHSkVSaOVQIqTclVDRWju/Pz6zYe3uzM7MNd/rmnk/H495mJ3r\n2rk+bXvtZ753ny+IiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIj4VDngdtdBiIiIv3zqOgAREUmN\ncjGe1xS4ANgG1ARqRx4/eBSXiIg4UibG82YB4SJeb5e8UERERESC51xgSpzfeyzwURJjEfFErF1J\nAOcDFwFnAm0jj/e9CErEkcrAMmAjkBN57SBgObAK+B/wCnAvsDaO918PXApsBb5KNFgR10YBzwPf\nAdnAYmCM04hEvHEOsAEbSwN4Gngt8vzPJP4H/QpgWoLvIeILBZ+evoj8Wxn40FEsIl4bB7wMhIBN\nwCGR1x8Ani10bj5wE/A11hIYABwJzAF+BiYC5fc4vw7wW6HXRAJpXuTfudgvdkWseS2SjqphM+42\nAt33eP1V4K+Fzs3HxhwqY7P3fgdmAllAFeBLoFuh79kCHJPsoEWSpWyM5/0bqA4MBj4DvgUmeBST\niGs/Y3/QD2DvgeZq2JTtwh4DfgGWYK3rN7F7ZGvk+QmFzt8WeS8RX4o1MTwGbAYmY5+EGgMDPYpJ\nxLWrgPrAu8Cje7y+GWsFFLZ+j+fbi/i6cqHzD8KSj4gvxZoYPt7j+Q7sl/rjYs4VCbJDgMeBnsCN\nwGVAm8ixL4BGCb5/HaACNvtJxJf2i3K8NnAYUAk4EVsQF8Y+NVXyNjQRJ/6JdR8VTMW+G/gXcBww\nHRtMjqZMMc/Bpnn/F9iZWJgi3omWGM4BrsE+5Qzd4/VtwH0exSTiykVAa2wQucAYbIrp3yOPLUBL\ndk/IKKoiQLjQ8z2/vhKbAisSeJfG+X1jsf7WnGKO3wksiDxygDw0KCf+djZa+SwCQF+s+6gM9gnq\nc6w0QDSnYTMyiksMezofG+wTEZEAKFjYVlAn5hjsU34ssogtMbwMXFfqyEREJKlinZVUMIDWEXgB\nK4mRTJWwpDM5ye8rIiKlFGti+Ax4G+gAvIV1K+UnMY4LsBIbmtstIuJYtFlJBa4Djge+weq81AB6\n7HG8GbZSNF5dKWEldffu3cNZWVl/fB0KhQiFQglcTkQk/ZUpUybWPXf2/r4kXX8B+y77L5CFVZNs\nXszxqsAKoC62SrQo4XC4qFmBIiJSnHgTQ6wthnhNwBb01ATWYCW7C6pKjor8exEwg+KTgoiIpFAq\nWgzJoBaDiEgpvPMOnHNOfC2GWAefRUQkIPLy4Oab4//+ZCWG35P0PiIikqCXXoI6deL//niaGUcB\nl2MziZrFf+lSUVeSiEgMdu6Exo1h7FgIhbztSqoD3AF8ii1uK4clBhER8ZEXXoD69aFt2/jfI1o2\nuQFrHRyCbYg+CXgDaBD/JeOiFoOISBQ7d8LRR8P48XDaad5NV/0nttL5NmBRPBcQEZHUGD8ejjjC\nkkIiYtmopwvwJLtbDeVL/A4REUm53Fx46CEbeE5UtDGGTdimIm2xTXu2YPsrLAUeTvzyIiKSDM89\nZ91Ip56a+HtF6386DPi+iNcbYYPPAxIPISYaYxARKcaOHZYUXnkFWrXa/bpXtZLeBA4G3sPGGj7E\ndllLNSUGEZFijBwJn38OY8bs/bqXRfQOAELAedh+uGuwhPEWsDqei8ZBiUFEpAi//QYNG8Ibb0CL\nFnsfS2V11SOAq4BWWHG8lvFcuJSUGEREijBkCMyZA5OL2OYsFYnhRGxNw2XASmy3tVFAbjwXLiUl\nBhGRQrZtg6OOgpkzoVkRdSi8WsdwNJYM/gJsxBa4lcG6lkRExKHhw+Hss4tOComIlk3ygX8Dt7J7\nPGElWvksIuLUTz9Bo0Ywd661GooSb4sh2jqGS7ANdD4AngHOJHl7OIiISJyGDIFLLik+KSQi1j/y\nlYELsW6ldsDzwBTg7eSHVCS1GEREItavh6ZNYcECOPzw4s9L5aykg4HO2AK3M+K5aByUGEREIvr2\nhXAYnnii5PNSmRhcUGIQEQFWr4Zrr4UXX4RatUo+16sxBhER8ZF//ANOPjl6UkhEtOmqIiLiE0uW\nwLRp8NVX3l5HLQYRkYC4/364+26oVs3b62iMQUQkAObOhS5drLVwwAGxfY9fxxjGYvs35JRwTghY\ngO0lPcvjeEREAicchnvugf79Y08KifA6MYwD2pdwvBowErgAOAabBisiInuYMcPWLnTvnprreZ0Y\nZgObSzh+BVaM77vI15s8jkdEJFDy8621MHAg7Jei6UKuB58bsnsjoPnA1W7DERHxl1degf33h4sv\nTt01XU9XLY+V8z4TqATMAeYCXxc+sX///n88D4VChEKhlAQoIuJKbq7NRBo9GuIbRo6P68SwBus+\n2s7uYn3HESUxiIhkgpdfhmOPhXbtUntd111JrwNtgHJYi+FkYInTiEREfGDrVhtbyM5O/bW9bjFM\nANpiW4CuAbKx7iOw3d+WYntHf4Ht/fAvlBhERBg8GM49F44/PvXX1gI3ERGfWbvWupAWLoR69eJ/\nH1VXFRFJEz17Qs2a8Mgjib2PV3s+i4hICi1enJpCeSVxPfgsIiJ76NcP7rsPqlZ1F4NaDCIiPjFz\nJixbBlOmuI1DLQYRER/Iz4c774RBg6BCBbexKDGIiPjAxImWEDr7oJSoZiWJiDj222/QtClMmACn\nnJK89/XrfgwiIhLFkCHQsmVyk0Ii1GIQEXGoYDHb/PnQoEFy31sL3EREAuiaa6B2bRt0TjYtcBMR\nCZj58213tmXLXEeyN40xiIg4EA7D7bfDgw9ClSquo9mbEoOIiAOvvQbbtkGPHq4j2ZfGGEREUmzH\nDmjSBMaMgTPO8O46mq4qIhIQw4fDccd5mxQSoRaDiEgKrV8PzZrBnDnQsKG319J0VRGRAOjVCw46\nCIYO9f5amq4qIuJzCxfCihU28OxnGmMQEUmBcBj69IEuXaBaNdfRlEyJQUQkBSZMgF9/tW07/U5j\nDCIiHvvlF2jcGF59FVq3Tt11NfgsIuJT99wD338Pzz+f2usqMYiI+NDXX1s57ZwcK5aXSn5d4DYW\nWA/kFHM8BGwBFkQe93scj4hISvXtC/36pT4pJMLr6arjgBFASQ2o94FOHschIpJy//43fPMNTJni\nOpLS8brFMBvYHOWcoHRniYjEbMcOay0MH257OQeJ6+mqYaA1sAiYDjR1G46ISHI8/riVvmjf3nUk\nped65fPnQD3gN+A8YCrQqKgT+/fv/8fzUChEKBTyPjoRkTisWQOzZsEzz7iOJD6p6MbJAqYBzWM4\ndyXQAvip0OualSQigXHppbaPc3a22ziCWivpUGAD1qXUEktUhZOCiEhgTJ8OixbBSy+5jiR+XieG\nCUBboCawBsgGykeOjQI6AzcBeVh3UleP4xER8cz27dC7N4wcCRUruo4mfkGZEaSuJBHxvexsWLIE\nJk1yHYnRymcREYcKVjgvXAh167qOxvh15bOISNoLh+HWW60mkl+SQiKUGEREEvTaa7B2Ldx2m+tI\nkkNdSSIiCdi2DZo0sf0WTjvNdTR70xiDiIgD998Pv/8Ogwe7jmRfSgwiIin22WfQoQN8+SXUrOk6\nmn1p8FlEJIXy8qBXL3jsMX8mhUQoMYiIxGHECKhWDbp1cx1J8qkrSUSklFavhhNPhDlzoGFD19EU\nT11JIiIpEA7DLbfYXgt+TgqJcF1ET0QkUCZPhhUr7N90pa4kEZEYbdkCTZvCK69Amzauo4lOXUki\nIh57/HHo1CkYSSER6koSEYnB++/DmDGQk+M6Eu+pxSAiEsVvv0HPnvDUU1C9uutovKcxBhGRKO66\ny/ZxnjjRdSSlE9StPUVEfO3TT+GFFzKjC6mAupJERIqRmwvXXWeDzn/6k+toUkeJQUSkGIMGQf36\ncPnlriNJLY0xiIgUYfFiaNcOFiwI7q5sWscgIpIkeXlw7bUwcGBwk0IilBhERAp54gmoXBmuv951\nJG5oVpKIyB6WLYOpU2HcOIivIyb4vG4xjAXWA9Emev0ZyAMu8TgeEZFi5eXBNddA165w1FGuo3HH\n68QwDmgf5ZxywKPAWwRnMFxE0tCQIVCpEtx0k+tI3PK6K2k2kBXlnN7Aa1irQUTEiZwcGDoU5s+H\nshk++ur6P78OcCHwdORrzUkVkZTLzYXu3eGRR2zdQqZzPfg8HLgHSwhlKKErqX///n88D4VChEIh\nj0MTkUwxcCDUrm1TVCU1ffpZwDSgeRHHVuwRQ03gN+B64I1C52mBm4h44rPPoEMHW8h22GGuo0mu\noBbRO2KP5+OwBFI4KYiIeGLHDujWDYYNS7+kkAivE8MEoC3WGlgDZAPlI8dGeXxtEZESPfAANGmS\nebWQognK9FB1JYlIUn34IXTpAl98kb6VU1UrSUQkRtu2Qd++8K9/pW9SSIRaDCKSca691spdjBnj\nOhJvBXXwWUQkpSZPhg8+sFlIUjS1GEQkY3z/PZxwArz+OrRq5Toa72mMQUSkBPn5ViDv5pszIykk\nQolBRDLCiBE26Py3v7mOxP/UlSQiaS8nB844A+bOhSOPdB1N6qgrSUSkCDt2QO/e8OijmZUUEqEW\ng4iktT59YPNmeP75zNuRTdNVRUQKeeMNeyxYkHlJIRFB+VGpxSAipfLdd9CiBUyZAq1bu47GDY0x\niIhE7NoFV11l3UiZmhQSocQgImln4EDbnvOee1xHEkwaYxCRtDJ7Njz1lG3AU66c62iCKTAthrw8\n1xGIiN/9+KN1IY0eDXXquI4muAKTGN5803UEIuJn4bCVvOjZE84/33U0wRaYxJDu5XFFJDFDh8KG\nDdCvn+tIgi8w01WrVg2zdCnUquU6FBHxm48/hosvhnnzoH5919H4R9pPV73kEnjhBddRiIjfbNoE\nXbtar4KSQnIEpsUwe3aY66+HJUu0glFETH4+dOwIzZvDY4+5jsZ/0r7FcOqp9kswZ47rSETELx59\n1EppDxzoOpL0EpjEUKaM7dOqQWgRAXjvPXj/fZg4EcqXdx1NevG6U2Ys0BHYADQv4viFwAAgP/K4\nC5hZxHnhcDjMunXQpAmsWQOVK3sVsoj43XffwZ//bBVTzz7bdTT+5deupHFA+xKOvwscB5wAXAM8\nW9Kb1aoFp58Or76atPhEJGByc6FLF9tjQUnBG14nhtnA5hKO/7rH88rApmhveO21MHZsomGJSFD9\n9a9wyCGqg+QlP4wxXAT8D3gT6BPt5A4dYPlyWLrU87hExGdefBHeegvGj7cieeINP/xopwJNgAuA\nqCsVypeHG2+E//s/z+MSER/54gu4/Xa796tVcx1NevNTddXZWDw1gB8LH+zfv/8fzxs1CnHbbSH6\n9oVKlVIWn4g48vPPtsh1+HBbsyDeSsVSsSxgGkXPSjoSWAGEgROBSZHXCttnB7cLLoALL7SCWSKS\nvnbtsnu9ZUt44AHX0QRLvLOSvE4ME4C2QE1gPZANFMw4HgXcDXQDdgK/AHcAnxbxPvskhrffhrvu\ngoULtRJaJJ3dey988gnMmKH1CqXl18SQLPskhvx8aNoUnn3WprCKSPqZMAHuuw8+/RRq1nQdTfD4\ndR2DZ8qWhVtvhREjXEciIl5YsMD2bJ46VUkh1QLbYgCrkVK/PixaBPXqOYhKRDyxYYONKQwebIvZ\nJD4Z12IAOOgg28bv6addRyIiyZKbC507w9VXKym4EugWA8BXX0GbNrB6NVSsmOKoRCTpbr7ZaiFN\nnapFbInKyBYDQKNG0KKFVVgUkWB76inYuNFWOCspuBP4FgPA9Olw//3w2WeauioSVG+/Dd26wUcf\nwZFFrWaSUsvYFgNA+/aQlQWzZ7uORETisWSJjRdOmqSk4AdpkRjKloVzz9XWfiJBtHGjVTIYPBhO\nO811NAJp0pUEsGMHNGhgqyOPPTZFUYlIQn7/Hc46yxLCww+7jib9ZNzK56I88ggsXmwDVyLib+Ew\n9Ohh65EmTdJgsxeUGIAtW+CII2wQOivL+6BEJH4PPmhjC6NHw4EHuo4mPWX04HOBqlXh+uth6FDX\nkYhIScaPh3HjYNgwJQU/SqsWA8APP0CzZrbD2yGHeByViJTaf/8LV1wBs2ZBkyauo0lvajFE1K5t\ny+hVXE/Ef3Jy4PLLbUxBScG/0q7FALYndChk/ZdVqngXlIjE7rvvoHVrm1betavraDKDWgx7OOoo\nOOcctRpE/GLrVujY0UrlKyn4X1q2GMCK6516Knz9tTYOF3Hp998tKbRpA9nZKluTSpquWoQePWyf\nhgEDPIhIRKLatWt3C2HiRChXzm08mUaJoQgrV8JJJ1nroUYND6ISkWKFw3DTTdZqnz4d9t/fdUSZ\nR2MMRWjQwGYoDR7sOhKRzNO/v+3VPHWqkkLQpHWLAWDNGjjuOPjf/+DQQ5MclYgU6Z//hCefhA8/\n1Hoil9SVVII+faxvc9iwJEYkIkWaOBHuvNPK4Ddo4DqazKbEUIKC1dA5OVCnThKjEpG9zJhhSeGl\nl1Tl2A/8OsYwFlgP5BRz/EpgEfAF8BHgya9S7dpw3XXw0ENevLuIAHzwAVx9NYwapaQQdF4nhnFA\n+xKOrwBOxxLCg8CzXgXSrx988gksW+bVFUQy17x50LkzTJhgq5sl2LxODLOBzSUcnwNsiTz/BKjr\nVSA1a1rhrjvv9OoKIplp0SLbgW3sWDjzTNfRSDL4abrqdcB0Ly/Qu7fNTnrnHS+vIpI5li6F886z\nWUjnn+86GkkWvySGdsC1QD8vL7L//ram4Y47IC/PyyuJpL+VK60m2aBBtl5I0kcqZiVlAdOA5sUc\nPxb4P2wsYnkx54Szs7P/+CIUChEKheIKJhyGM86Av/wFbrwxrrcQyXirVsFVV9njhhtcRyPF8fN0\n1SyKTwyHAzOBq4C5JbxHQtNVC1u4ENq3t2awCuyJlM6qVdCuHdx+u3XPin/5NTFMANoCNbFpq9lA\n+cixUcBo4GJgdeS1nUDLIt4nqYkBoGdPqF5d5TJESkNJIVj8mhiSJemJYd06OOYYmDvX9m8QkZKt\nXm0bYPXta9UExP/8usDNt2rVsqmrajGIRLd6tbUUbrtNSSETZGxiAGsOz5oFU6a4jkTEvwq6j/r0\nscQg6S9ju5IKvP8+XHml9ocWKcry5XDWWXD//TYuJ8GiMYYE9OwJFSvaIh0RMV9+CeeeCw88AL16\nuY5G4qHEkIDNm6366uTJcMopnl1GJDA+/9z2aR482NYqSDBp8DkB1avbXg29ekFurutoRNyaM8fW\n+YwcqaSQqZQYIi67DA4/HIYMcR2JiDszZ0KnTjB+PFxyietoxBV1Je1h1Spo0cI+MTVs6PnlRHzl\n9ddtO86//93WK0jwaYwhSYYNs0VvEyZAWbWnJEOMHQt/+xtMmwYnneQ6GkkWJYYkycuzmvIXX2wr\nPEXSWTgMjz0Gzzxj23I2auQ6IkkmJYYkWrECTj7Z+lubF1cTViTg8vPhrrssIcyYof3Q05FmJSXR\nEUfYNL0rroAdO1xHI5J8ublw993WbfrBB0oKsje1GIq9oM1UqlcPHn88pZcW8dTPP9uMoxo14Lnn\n4MADXUckXlGLIcnKlLF+10mTtBWopI9vv4VTT4Vjj4WJE5UUpGhKDCWoUQPGjYMePeDHH11HI5KY\n+fMtKfTqBcOHQ7lyriMSv1JXUgzuuMPGGkaOtJaESNBMmwbXXQfPPgsXXeQ6GkkVzUry0I4dcPrp\ntk/0X//qLAyRuDz9NDz4IEydCi2L2h9R0la8iWG/ZAeSjipWtLGGk0+2ldFaFSpBkJtr+yd8+y18\n9BE0aOA6IgkKtRhK4d13oVs3mDcP6tZ1HY1I8TZuhM6doWpVePFF7TWSqTQrKQXOOss2QO/SRVVY\nxb8WLbIuozZtrPtISUFKSy2GUsrPtzngdetqYx/xn8mT4cYbYcQI6NrVdTTimgafU2jLFhtvyM6G\nyy93HY2IfWAZMMCmV0+ZAiee6Doi8QMNPqdQ1apWovj006FWLdsoXcSVzZvhmmvgoINs/OvQQ11H\nJEHn9RjDWGA9kFPM8cbAHGAHEKiJoEcfbaW5u3aFJUtcRyOZasECK5OdlWWls5UUJBm8TgzjgPYl\nHP8R6A0Ect+0M86wYnsdO8K6da6jkUwzZgyccw48/DA88QRUqOA6IkkXXnclzQaySji+MfLo6HEc\nnunWDVauhAsugFmzVHtGvLd9O9x6q+00+MEH0KSJ64gk3Wi6ahI88AA0a2bdSnl5rqORdLZkCbSP\ntMHnzVNSEG8oMSRBmTJWg+aww+D6622GiEgyhcP2O9a2rbVSR4+GypVdRyXpKjCzkvr37//H81Ao\nRMhndSkqVLB9G9q3h5tvtvo0KrgnyfDTT/aB45tvYPZsaNzYdUSS7lLxpysLmAaUtElmf2AbMLSY\n475ax1CSrVttQLBVKxg2TMlBEjN7Nlx1le1B/sgjVrdLJFZ+XeA2AWgL1MSmrWYD5SPHRgG1gE+B\nKkA+lhyaAr8Uep/AJAawHbLOPNNKaDzyiJKDlF5uLjz5JAwdat1GHQM7PUNc8mtiSJZAJQawjX3a\ntbPyGXv0golElZMD3bvb4skxY6B2bdcRSVCpiJ7P1Khh1Vg//dRmLQUsr4kDeXkwaJCtj7nlFvjP\nf5QUxA21GDy2YYMNSJ9yihU2K6tULEVYtsxaCQceaCuY69d3HZGkA7UYfOqQQ+C992DxYhtE3LnT\ndUTiJ3l5NoPt1FPh6qvhnXeUFMQ9tRhSZPt22xo0Lw9eew0qVXIdkbg2bx7ccIN1O44aBUce6Toi\nSTdqMfjcAQdYrfyjjrLyGaqtlLm2bLENnzp1sj3E33lHSUH8RYkhhcqXt2JnoZDt57BokeuIJJXC\nYWstNmtmLcgvv7TuRU1nFr8Jyq9k4LuSCnv1VZt5MmaMfXKU9LZypbUSVqywbqPTTnMdkWQCdSUF\nzGWX2XTEv//dpiiqvlJ6+uUX+38cCtkA88KFSgrif0oMDrVsCW++aY8OHWDTJtcRSbLk58Pzz1td\noxUr4MMP4d57tWeCBIMSg2OHHQYzZ8Lxx9s+vR9/7DoiSdTHH1utrJEjYdIkeOklqFfPdVQisdMY\ng49MmwY9e0J2tk1jLFfOdURSGkuXWrfR6tW2kc6VV2pBo7ilMYY0cMEF8MknMGWK9Ul/843riCQW\na9ZYQm/TBlq0sBbg1VcrKUhw6VfXZ7Ky4K23rMzyySfDU09pYNqvNmyA++6D446DP/0Jvv4a7rlH\n27tK8Ckx+FC5cnDHHTZgOX48XHqpDWCKP6xdC7ffbgPLFSrYeoRBg6B6ddeRiSSHEoOPNW4MH31k\nG/+0bAn/+Afs2OE6qsy1apXtzte8uS1KW7zYSqqrAqqkGyUGn9tvP7jpJvj8c6vTf8wxtv5BUmfR\nIujXz2aNVa1qg8yPP24zykTSkWYlBcxbb9mOXj//bLvDnXSS64jS065dloCHD7eS2H362L7LBx/s\nOjKR2GkHtwyyc6eV0hgwAE4/HR56yIrzSeK2boXnnrNtNatXt7GEzp21ME2CSYkhA/36q32i/e9/\n4fDD4e67oWlT11EFTzhs04RHj4YZM6x0Re/e0Lq1CtxJsCkxZLDNm21a64gRNkjdr5/9cZOSbdxo\nq5JHj7ZB/Z49bRc1DSZLulBiELZvt26QIUNsRlO3bnDRRbD//q4j84+tW+H112HiRJt2eswxlhDa\ntlXrQNKPEoP8IS8Ppk+3vR9ycmwVbs+e0KSJ68jc+PVXG0h+5RV4911LAl272krzgw5yHZ2Id5QY\npEjLl9tA9XPPWWG3tm1tMLVuXdeReWvVKksG//kPfP+9rUzu2tVWlGshmmQKvyaGsUBHYAPQvJhz\nngTOA34DrgEWFHGOEkOCdu60T8uTJllXytFHW0siFLJup6B3o2zbZosBP/8cXn4Z1q+H886Djh3h\n3HOhWjXXEYqknl+L6I0D2pdwvANwFNAQ6AU87XE8+5g1a1aqL5mQeOMtX97+UI4dCz/8YFVA166F\n9u2tJHSPHvYH1Yu9qL34Ga9fb91ld91lA+61a9u6jooVbTB53TrbD+Evf4kvKWTK74VLijklQvF8\nk9eJYTawuYTjnYDxkeefANWAQz2OaS9B+x+djHgrVLAk8dBD8O23Vg30pJNsu9Hu3aFBA+jSxf7Q\nvvuuzXpyFfOuXbYt5rRpVn6iUyeoU8fGS15/HapUscH2TZtg1iyrMdWqVeIlyzPx9yLVFHNKhOL5\npv2SHERp1QHW7PH1d0BdYL2bcDJPmTLQqJE9brnFKrkuXw7z59tjwABYsMDm9OfmwhFHWOIo+Lde\nPahRwz6px9No3bnTPv2vW2ctmZ9+shpEX31l1UpXrrT3P+kkW6PRvbsNqmdlBb/7S8SvXCcG2Hec\nQ4MJDpUtuztRXHGFvbZrl7UsVq60Kq8rV8Ibb9jz6tXtkzrY82rVLFns3GnvVbasnbdsmQ0C//LL\n3o8qVez9a9Wy7qATTrBE0K2bxXDkkVCpkqufhoh4JQvIKebYM0DXPb5eStFdSQuxhKGHHnrooUfs\nj+fwqSyKTwwdgOmR562AuakISERE3JkAfA/kYmMJ1wI3RB4F/gksBxYBJ6Y6QBEREREREdpj4xdf\nA/2KOSeELa5bDMxKSVQlixZzTeAtbPxlMbY40KWx2Gyy4roNwRYzfo21Dk9IRVBRRIv5SizWL4CP\ngGNTFFdxYvkZA/wZyAMu8Tyi6GKJOYR/7r1o8frtvgOoB7wHfInF1KeY8/x2/zlVDuuqygLKY/9D\nC1cNqob9UAsKRdRMVXDFiCXm/sCgyPOawI+4nWV2GvbLFst40sn4YzwpWsynAFUjz9vjPuZo8YL9\n7swE/g1cmoqgoogWs9/uvWjx9sdf9x1ALeD4yPPKwDL2/XtRqvsvE7b2bIn9kf0W2AlMBC4sdM4V\nwGRsHQWi/yKDAAACy0lEQVTAplQFV4xYYv4BqBJ5XgX7Bc1LUXxF8f1ixiJEi3kOsCXy/BN2//Fy\nJVq8AL2B14CN3ocTk2gx++3eixav3+47gHXYh0eAX4D/AYU3ni3V/ZcJiaGoRXR1Cp3TEDgYa47N\nB65OTWjFiiXmfwHNsMH9RcBtqQktbsUtZgyK69j9icuv6mAfIApKy4QdxhIrv9170fj9vsvCWjyf\nFHq9VPef6yZQKsRyc5THZkSdCVTCPinOxfrjXIgl5vuwTwkh4EjgHeA4YJt3YSUsqIsZ22Ez6vy+\n/dFw4B7s51qGYFRP9tu9F42f77vKWGvxNqzlUFjM918mtBjWYoMzBeqxu9laYA3wNrAdaxp+gP3P\ndiWWmFsDkyLPvwFWAkd7H1rcCv831Y285nfHYp8SOxG9G8e1Fli340psfOEpLG4/89u9F41f77vy\nWJfci8DUIo4H9f7zzH7Y/8AsoAJFD+Q2Bt7FBu4qYQNPLndPjiXmx4HsyPNDscRxcIriK04WwVvM\nmEXxMR+OjfW0Slk00WURfVYSWGVjP8xKgpJj9tu9ByXH68f7rgzwPDCshHP8ev85dR42Ur8cuDfy\nWuGFdndisyNyKH66VypFi7kmMA3r58zBBvFcCuJixmgxj8Y+xS6IPOY5iHFPsfyMC/glMcQSs5/u\nvWjx+u2+A2gD5GMfIAt+V8/D//efiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIgEXz7wwh5f74cVVZwW\n6xtkQkkMEZFM8itW6K9i5OuzsRXaMdcmU2IQEUk/04GOkeeXYyu6Yy6qqMQgIpJ+XgG6AvsDzdm3\nDHeJlBhERNJPDlYM8HLgP6X95kzYj0FEJBO9AQwB2gJ/Ks03KjGIiKSnsdgeIl9iGwvFTF1JIiLp\npWD20Vqs1HbBa0HZMVFERERERERERERERERERERERERERERERERERDLV/wO/Iw1z0TEqvQAAAABJ\nRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x86a5db0>"
       ]
      }
     ],
     "prompt_number": 24
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Safe to say they are they same!\n",
      "\n",
      "-------------"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "References\n",
      "----------\n",
      "\n",
      "[1] : http://www.potto.org/gasDynamics/node80.php"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Question 3\n",
      "=========="
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "T0,P0,P,A,Pamb,gamma,M,V = sp.symbols('T_0 P_0 P A P_amb gamma M V',real = True,positive = True)\n",
      "P0_eq_temp = sp.Eq(P0/P,(1+(gamma-1)/2*M**2)**(gamma/(gamma-1)))\n",
      "Math(sp.latex(P0_eq_temp) + r\"~~~~~\\text{where}~~ \\gamma = 1.667\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "$$\\frac{P_{0}}{P} = \\left(M^{2} \\left(\\frac{\\gamma}{2} - \\frac{1}{2}\\right) + 1\\right)^{\\frac{\\gamma}{\\gamma - 1}}~~~~~\\text{where}~~ \\gamma = 1.667$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 60,
       "text": [
        "<IPython.core.display.Math at 0x86b6670>"
       ]
      }
     ],
     "prompt_number": 60
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The only unknown is $M$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Pamb = 1*1013*1000\n",
      "vals = {T0 : 500, P0 : 2*1013*1000, A : 25.0/100/100, P : Pamb, gamma : 1.667}\n",
      "P0_eq = P0_eq_temp.args[0]-P0_eq_temp.args[1]\n",
      "P0_eq.subs(vals)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "$$- \\left(0.3335 M^{2} + 1\\right)^{2.49925037481259} + 2$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAUgAAAAgBAMAAACC6vyDAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEM3dMkS7mXYiie9U\nZqsqREkJAAAFKElEQVRYCeVXa2gcVRT+ZndmXzOTnbbQav2x21WwStEVX41VWQStqJD94QMF7VZ8\nIApOW0lqpXUUoRrBzq8UhbiL/mgaH41iS3/42IKPCoGsIGoI2K1UUdKSjmmokTTruXd2dmaSbLLJ\nbvPHA/fMved859xv75x77yzQcol1buM5Q8B3/XY78LQmnuotobsbiHWfxIcrVhAkB7WPoGIOcAAs\nYFnkIH7j8wxpkXJEZy1mtiWlyjTUglLGbnRirFLJQbHwI1ZD3JOGA2DgZeGIdnSYNJNwSItnoxZr\nQVOclL7WEMzAwi6cwO8PPKhhrYVnEACUNBwAAy8PycPIZ2mmSELLZzDBmpyNTkhkS5j4nDzbUcQ3\nwCsWdmTe5SQdAAO3kCQVUf1sQ8x3X0JLZDDFGhGxpJ17tQ4Nt9ESf0atAPGYhVXnC5ykA7DB9RMv\nzsOKqH7ERuYqJDS5KF5gDQgagqnkiORGCK8aQJuGmGBBvba9StIG2OD6ibnnzgX8NTcrotrA7Ri8\nG9PpoZoJTclFJlij3UsWYZJq9Sh1OoG7gGEi+ZG4SWM16QBsMA3nEXYYNCasiOZAPsZtPzN9CRIa\n7n5tkjdVR8gUJuhtsqiOIvohGIIlJBEocpJVAAez6HkkME+dzQyjA2S2xEpkkwpSho7Hrpd3C5By\nrOER3B80aZdnYVFF5ktYD6Wr6/xwDhGbZBXAwbPT+iw7faP5B0FjLv91ZLxj8Futh35uXgsZoQxr\n4obBrTEtpkd0OifbMaLhC8LRaXMVK05lAA6AgedK67GxM2pN3wvc8k7/c4Cr1mubs1RnY7RRgX1H\ninYRcWRNraPeXmpfVSq414T6/Gmp/0WwJlcqZ9G/i4IOaLin+y/gRsL9NF1c03cQ0R3ThgNg4Bmy\nv/MJn4UV/AasMskY1TFUchVN/A+Dtv1LShxjfAnrlx72/vN+W/MjoYxDNJ0rchHqAGJJskSmEE+7\nClf/yaij90tS6hgpKiJmoD3LHcDqPYxksMSNLVCX2jmoHOQBbzo6xkJldqVS8R9GIukqlG3c8OP0\n/OEIrSZVmW1SDPuJCIsL56qjph/H7QxyGm3nvMlGgHgZEq870Osmn6PKNk5/iSrhYXrzvIi4zU8y\n4vvVdtDSdJVkwJpBcjvdrjokXn3AlSy5o878wveTvo8qQfX9Mj9JiS1nS6RKknIFfTmfosKnO5Vt\nDojv6161lX80REtDGRxrO+tl4Scp+n6AF7fYvktyS8Ebe42HJPBBhnyuCqVpy5j5olCUc94gP8k5\n7yEvvOG+S/ImYOX1TK5gwTfT6y7UXrf8CZlcFaa7Ooy4EUHeYGBHZpCsFrTjXfqzRnLGXqSVjOtQ\n2DxRE6EpV9EnT4xu3TBk/XucyNSmjqVSl3+aSlFlUK2y0mnRSr6XSm1KpfiGwFu12XiHSAZzUNlk\ncQuhC66SkwiTuReBdSWMcrCj/Cs557eRA13U01lJRcfb3sA/aDkGEE6STc4hOOmqMNEn8zDCRwE6\nJj3iJ9n63b0ZeMMzHd6kwa24LCP8DbWEIcNVdEGOmEAakXMQq0dUNdJPkm6sFslxO0/0ln76o+aR\nLRqw9vVTwO3AQ11PetX+j7cBo+OGYijPjp/0xMAhGR4dP3MRbpxApVLxkZSz3tkb7DskbXjAaDBs\nQVh1JWfhwoVZpoUNUsmLoa/xFsmjdfKIA3UcjZt/bRy6VOTppQbW4m6o9S5aZ2WzmdVCsxkWjqeD\nsznpaS68sejBxmB1UXRO/W/lP0XhoQFBEnw0AAAAAElFTkSuQmCC\n",
       "prompt_number": 52,
       "text": [
        "                 2.49925037481259    \n",
        "  \u239b        2    \u239e                    \n",
        "- \u239d0.3335\u22c5M  + 1\u23a0                 + 2"
       ]
      }
     ],
     "prompt_number": 52
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "M_a = float(sp.nsolve(P0_eq.subs(vals),2))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 27
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Math(\"M = \" + sp.latex(M_a))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "$$M = 0.978965655189$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 28,
       "text": [
        "<IPython.core.display.Math at 0x8624130>"
       ]
      }
     ],
     "prompt_number": 28
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "T0_eq_temp = sp.Eq(T0/T,(1+(gamma-1)/2*M**2))\n",
      "Math(sp.latex(T0_eq_temp) + \"~~~~~\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "$$\\frac{T_{0}}{T} = M^{2} \\left(\\frac{\\gamma}{2} - \\frac{1}{2}\\right) + 1~~~~~$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 77,
       "text": [
        "<IPython.core.display.Math at 0x3832b50>"
       ]
      }
     ],
     "prompt_number": 77
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Solve for T"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "T0_eq = T0/T - (1+(gamma-1)/2*M**2)\n",
      "T0_eq = T0_eq.subs(vals).subs(M,M_a)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 30
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "T_a = sp.solve(T0_eq,T)[0]\n",
      "Math(\"T = \" + sp.latex(T_a) + '[K]')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "$$T = 378.897630800344[K]$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 74,
       "text": [
        "<IPython.core.display.Math at 0x86dc4d0>"
       ]
      }
     ],
     "prompt_number": 74
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Rhe = 2.08*1000; R = sp.symbols('R') \n",
      "V_eq = M*sp.sqrt(gamma*R*T)\n",
      "V_a = V_eq.subs({M:M_a, gamma:1.667,R:Rhe, T:T_a})\n",
      "print 'a)'\n",
      "Math(\"V = M\\cdot c = M\\cdot \\sqrt{\\gamma R T} = \" + sp.latex(V_a)+'[m/s]')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "a)\n"
       ]
      },
      {
       "latex": [
        "$$V = M\\cdot c = M\\cdot \\sqrt{\\gamma R T} = 1122.09045495887[m/s]$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 75,
       "text": [
        "<IPython.core.display.Math at 0x7d09650>"
       ]
      }
     ],
     "prompt_number": 75
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "m_dot_a = V_a*0.0025*1013000/(378.9*2080)\n",
      "# m_dot\n",
      "print 'b)'\n",
      "Math(r\"\\dot{m} = \\rho\\cdot V\\cdot A = V\\cdot A \\cdot \\frac{P}{R\\cdot T} = \" + sp.latex(m_dot_a) + '[kg/s]')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "b)\n"
       ]
      },
      {
       "latex": [
        "$$\\dot{m} = \\rho\\cdot V\\cdot A = V\\cdot A \\cdot \\frac{P}{R\\cdot T} = 3.60569827281319[kg/s]$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 76,
       "text": [
        "<IPython.core.display.Math at 0x8614150>"
       ]
      }
     ],
     "prompt_number": 76
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "-----------------\n",
      "\n",
      "Maximum mass flow rate is achieved at $M = 1$. The area is fixed. \n",
      "\n",
      "Thus, we must find the ambient pressure which would allow sonic flow at the nozzle outlet. "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "P0_eq_temp\n",
      "P_c = sp.solve(P0_eq.subs({M:1,P0:2*1013*1000,gamma:1.667}),P)[0]\n",
      "Math(sp.latex(P0_eq_temp)+'~~\\\\Rightarrow~~P = '+ str(P_c/(1013*1000))+'[atm]')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "$$\\frac{P_{0}}{P} = \\left(M^{2} \\left(\\frac{\\gamma}{2} - \\frac{1}{2}\\right) + 1\\right)^{\\frac{\\gamma}{\\gamma - 1}}~~\\Rightarrow~~P = 0.974184338881403[atm]$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 91,
       "text": [
        "<IPython.core.display.Math at 0x86cbe10>"
       ]
      }
     ],
     "prompt_number": 91
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "T_P_rel = sp.Eq(T0/T, (P0/P)**((gamma-1)/gamma))\n",
      "# T_P_expr = T_P_rel.args[0]-T_P_rel.args[1]\n",
      "# sp.nsolve(T_P_expr.subs({T0:500,P0:2,P:P_c,gamma:1.667}),1)\n",
      "# sp.solve(T_P_expr,T)[0].subs({T0:500,P0:2,P:P_c,gamma:1.667}).evalf()\n",
      "Math(sp.latex(T_P_rel)+'~~\\\\Rightarrow~~T = '+ str(374.9247407)+'[K]')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "$$\\frac{T_{0}}{T} = \\left(\\frac{P_{0}}{P}\\right)^{\\frac{1}{\\gamma} \\left(\\gamma - 1\\right)}~~\\Rightarrow~~T = 374.9247407[K]$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 105,
       "text": [
        "<IPython.core.display.Math at 0x7d09230>"
       ]
      }
     ],
     "prompt_number": 105
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "V_c = V_eq.subs({M:1, gamma:1.667,R:Rhe, T:374.9247407})\n",
      "print 'c)'\n",
      "Math(\"V = M\\cdot \\sqrt{\\gamma R T} = \" + sp.latex(V_c)+'[m/s]')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "c)\n"
       ]
      },
      {
       "latex": [
        "$$V = M\\cdot \\sqrt{\\gamma R T} = 1140.17500802006[m/s]$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 109,
       "text": [
        "<IPython.core.display.Math at 0x86c2cb0>"
       ]
      }
     ],
     "prompt_number": 109
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "m_dot_c = V_c*0.0025*P_c/(374.9247*2080)\n",
      "# m_dot\n",
      "print 'c)'\n",
      "Math(r\"\\dot{m} = V\\cdot A \\cdot \\frac{P}{R\\cdot T} = \" + sp.latex(m_dot_c) + '[kg/s]')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "c)\n"
       ]
      },
      {
       "latex": [
        "$$\\dot{m} = V\\cdot A \\cdot \\frac{P}{R\\cdot T} = 3.60707130206581[kg/s]$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 114,
       "text": [
        "<IPython.core.display.Math at 0x87fc6f0>"
       ]
      }
     ],
     "prompt_number": 114
    }
   ],
   "metadata": {}
  }
 ]
}