MAE6286-Module1.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:6aeec53942fc13d2aea458ac24092a58dbd19c5ebedfee284fae6ff27ee6f55f"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "a = numpy.linspace(4,23,43)\n",
      "a[5]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 5,
       "text": [
        "6.2619047619047619"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "p = 7\n",
      "r = numpy.array([11.2, 4.7, 6.6])\n",
      "p/r"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 8,
       "text": [
        "array([ 0.625     ,  1.4893617 ,  1.06060606])"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "spacings = numpy.array([0.04, 0.02, 0.01])\n",
      "results = numpy.array([1.600, 1.500, 1.475])\n",
      "def convergence_order(spacings, results):\n",
      "    if spacings[2]/spacings[1] == spacings[1]/spacings[0]:\n",
      "        return numpy.log((results[2] - results[1])/(results[1] - results[0])) / numpy.log(spacings[1]/spacings[0])\n",
      "    else:\n",
      "        return 0\n",
      "convergence_order(spacings, results)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 10,
       "text": [
        "2.0000000000000062"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Rocket flight coding assignment\n",
      "%pylab inline\n",
      "import numpy\n",
      "import math\n",
      "import matplotlib.pyplot as plt\n",
      "\n",
      "# Given constants\n",
      "ms = 50. # Rocket shell mass, kg\n",
      "g = 9.81 # Acceleration of gravity, m/s**2\n",
      "rho = 1.091 # Average air densigy, kg/m**3\n",
      "area = numpy.pi*(0.5)**2 # Max cross sectional area, m**2\n",
      "ve = 325. # Exhaust speed, m/s\n",
      "C_D = 0.15 # Drag coefficient\n",
      "\n",
      "\n",
      "mp_dot0 = 20. # Propellant burn rate, kg/s\n",
      "tpmax = 5. # Maximum propellant burn time at mp_dot rate, s\n",
      "\n",
      "# Initial Conditions\n",
      "h0 = 0. # Assume that the rocket starts at height = 0, m\n",
      "v0 = 0. # Assume that the rocket starts at rest, velocity = 0, m/s\n",
      "mp0 = 100. # Initial propellant mass, kg\n",
      "\n",
      "# State vector is as follows:\n",
      "# x[0] - Height, m\n",
      "# x[1] - Velocity, m/s\n",
      "# x[2] - Mass of remaining propellant, kg\n",
      "x0 = numpy.array([h0, v0, mp0])\n",
      "\n",
      "# Equations of motion\n",
      "# Propellant burn rate is a function of time so the derivative of the state vector takes time as a parameter.\n",
      "def xdot(x, t):\n",
      "    h = x[0]\n",
      "    v = x[1]\n",
      "    mp = x[2]\n",
      "    hdot = v if h >= 0 else 0\n",
      "    mp_dot = -1*mp_dot0 if t <= tpmax else 0\n",
      "    vdot = (-1*(ms + mp)*g - mp_dot*ve - 0.5*rho*v*abs(v)*area*C_D)/(ms+mp) if h >= 0 else 0\n",
      "    return numpy.array([hdot, vdot, mp_dot])\n",
      "\n",
      "# Euler method\n",
      "def euler(x0, dt, tmax, xdot):\n",
      "    n = math.ceil(tmax/dt) + 1\n",
      "    t = numpy.arange(0, n*dt, dt)\n",
      "    results = numpy.zeros((len(x0), n))\n",
      "    results[:,0] = x0\n",
      "    for i in range(1, n):\n",
      "        results[:,i] = results[:,i-1] + dt*xdot(results[:,i-1], t[i-1])\n",
      "    return t, results\n",
      "\n",
      "dt = 0.001 # Time step, s\n",
      "tmax = 50. # Max time, s\n",
      "t, results = euler(x0, dt, tmax, xdot)\n",
      "\n",
      "plt.figure()\n",
      "plt.plot(t, results[0,:])\n",
      "plt.figure()\n",
      "plt.plot(t, results[1,:])\n",
      "plt.figure()\n",
      "plt.plot(t, results[2,:])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Populating the interactive namespace from numpy and matplotlib\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 52,
       "text": [
        "[<matplotlib.lines.Line2D at 0xb0496b0c>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEACAYAAAC6d6FnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XucjnX+x/HXREolhw6M05IU0wGpdO5WCFsOj98unXaF\n3Q4q2gqRjV8nVlLajd0Uq0JpK1HSCHfpwFQUYTDWOVS/0IFymPn98bmm+26aYWbuw/e67+v9fDzu\nx33f1336dDW+n+t7BhERERERERERERERERERERERERERSWMTgO3AsmJeuwvIB2pEHRsErAFygXZR\nx1t637EGGJOQSEVEJK4uBlrw6wRQD5gNrCOSALKAT4HDgQZAHpDhvZYDnOs9ngW0T1jEIiJyUIeV\n8n0LgB3FHB8NDChyrDMwFdgHrMcSQCsgE6iCJQGAZ4EuZQtXRETipbQJoDidgc3A0iLHa3vHC20G\n6hRzfIt3XEREHKhYzs8dBQwG2kYdyyjhvSIi4kPlTQCNsPb9z7zndYFPsKaeLVjfAFGvbfaO1y1y\nfEuxX96oUcHatWvLGZqISGCtBU5OxBc3oPhRQFB8J3AloKEXUGHtYBGWJDI4eCdwgZihQ4e6DsE3\ndC4idC4idC4igIKyFOql7QOYCnwAnAJsAnoWLbCjHq8Apnn3bwJ9ol7vAzyNDQPNw0YQiYiIA6Vt\nArrmEK+fVOT5w96tqE+AM0r5myIikkCxjAKSJAiFQq5D8A2diwidiwidi/Lz68gdrzlLRERKKyMj\nA8pQrpd3FJCkub17YckSWLQIVq+GNWtg0ybYudNuBw5ApUpwxBFwwglQpw785jdw5pnQooXdqlRx\n/V8hIgejGoD8bP16eO01u+XkwMknw3nnQdOm0Lgx1K8P1atD1apQsaIliR9/hC+/hC1bYN06WLrU\nEsfnn0Pz5tC2LXTubIkhw69/bSJpoqw1AL/+k1QCSJIffoBp0+Cpp2DtWrjqKiuwW7eO7Qp+925Y\nsACys+E//7Hvuu466NkTatWKX/wiEqEEIKXy5ZcwerQV/BdeCDfeCB062JV9vOXnwwcfwLPPwksv\nQadO0K8fnHVW/H9LJMjKmgA0CihgduyAu++GJk3g229h8WKYOdOu/BNR+AMcdhhcdFGklnHaaVbL\nuOoq+OyzQ39eRBJDCSAg9u+HceOs4P/hB1i2DMaOhQYNkhtHjRowYIB1KrdtC+3bw9VXw4YNyY1D\nRJQAAmHpUjjnHGvrnzPHEkEdx+uwHnkk9O0LeXmQlQUtW8LDD8NPP7mNSyRIlADS2P798NBDcPnl\nVtjOm2ejcfzk6KPhvvvgo49s5NGZZ8L777uOSiQY1AmcprZsge7doXJlmDAB6tU79Gf8YPp0uOUW\n+OMf4f77bZ6BiJSOOoGF+fOtyadjR3jrrdQp/AG6dLGO4TVr4OyzbT6BiCSGEkCa+fvf4dprbcjl\n4ME2AifVnHgivPwy3HmnzUd4/nnXEYmkJzUBpYn8fOjfH2bNgjffTP7onkRZuhR+9zu47DJ4/HHr\nPBaR4qkJKIB++smGUn78sU24SpfCH6xT+OOP4euvrTawfbvriETShxJAivvxR2s3z8+3ZReqV3cd\nUfwde6zNIG7XztYmWr7cdUQi6UFNQClszx4r/KtXt3byRM3k9ZPnn7e+geeegyuucB2NiL9oLaCA\n+PFHW0rhxBNh0qRgFP6F3nvP+gVGjYLrr3cdjYh/KAEEwIED0K2bFfpTpkCFCq4jSr4VK6wGMHAg\n3Hab62hE/CFRncATgO3AsqhjjwArgc+AV4CqUa8NwjZ+zwXaRR1v6X3HGmBMaYOUiIICuPVW2LXL\nhnoGsfAHWz5iwQIYMwYeeMDOi4iUTWkTwESgfZFj2cBpQDNgNVboA2QB3b379sBYIhlpHNAbaOzd\nin6nHML999uyCa+8olmyDRpYEnjpJVtgTklApGxKmwAWADuKHJsD5HuPFwF1vcedganAPmA9kAe0\nAjKBKkCO975ngS7lCTqoXngBJk60sf7HHus6Gn+oVQvCYZg715qDlARESi9ew0B7AbO8x7WBzVGv\nbQbqFHN8i3dcSmHxYrj9dtuusWZN19H4S40atsppdrbNflYSECmdeIwduRfYC0yJw3dJMbZvh65d\nbRnnZs1cR+NPxx0Hb79tM4YrVrSmMu1BLHJwsSaAG4COwOVRx7YA0cuP1cWu/LcQaSYqPL6lpC8e\nNmzYz49DoRChUCjGUFPT/v024qdHDxv6KCU7/nhrCmrd2paMuPde1xGJJFY4HCYcDpf782W5RmoA\nzATO8J63Bx4FLgW+jnpfFlYbOBdr4nkbOBkowPoK+mL9AG8ATwCzi/ktDQP1DB1qyzu89VZqLuzm\nwtattgXlgAFw002uoxFJnrIOAy1tDWAqVtAfD2wChmKjfiphncEAHwJ9gBXANO9+v3essDTvA/wb\nqIz1GRRX+Itn/nwYP97a/1X4l15mpiXMSy6xpiHVnESK59dW0sDXAL78Es46y0b9tG3rOprUtGSJ\nTRZ74QXrGxBJd5oJnAYKCmyNn6ZNYcQI19GktnDY+lBmz7aEKpLOtBx0Gpg8Gdats5EsEptQCP75\nT+jUCTZtch2NiL+oBuAzX3wBzZvrijXeRo2yFUTfew+qVHEdjUhiqAkohRUU2JVqixa6+o+3ggIb\nEbRli02mC9LqqRIcagJKYVOmwMaNMGSI60jST0YGPPkk7N1r+wmIiGoAvrFzp3X6Tp8OrVq5jiZ9\n7dwJF14IN99sS2uIpBM1AaWo226DffvgX/9yHUn6W7cOzj/falwaHirpRAkgBX3yCfz2t7bX7XHH\nuY4mGObOheuug4ULbVlpkXSgPoAUk58PffrA8OEq/JPp8stt+eiuXWH3btfRiLihBODYc8/ZMg89\neriOJHjuuANOPx3+9CctIS3BpCYgh3bvhlNPhWnTrE1akm/PHls47tpr4a67XEcjEptELQYnCfDY\nY1bwq/B3p3JlePVVG3nVvLk1DYkEhWoAjmzfbhub5+RAo0auo5F58+D66+Hjj6F2bdfRiJSPRgGl\niD59bFP3xx5zHYkUeuAB21Vs7lzNFJbUpASQAtautSaHVas08sdPDhyADh3g7LPh4YddRyNSdkoA\nKeCGG2zsedSul+IThfswPPUUdOzoOhqRslEC8Lk1a6zTNy8PqlVzHY0UZ8EC20Xso4+gfn3X0YiU\nniaC+dwDD0Dfvir8/ezii21IaPfutnicSLpSDSCJVq2yMed5eVC1quto5GDy821p7qwsGDnSdTQi\npZOoGsAEYDuwLOpYDWxD+NVANhB9TTsIWAPkAu2ijrf0vmMNMKa0QaaLBx+Efv1U+KeCww6z/Zgn\nT7aRQSLpqLQJYCLQvsixe7AEcAow13sOkAV09+7bA2OJZKRxQG+gsXcr+p1pa8MGmDVLSxCnkhNO\ngEmTbJmOr75yHY1I/JU2ASwAdhQ51gmY5D2eBHTxHncGpgL7gPVAHtAKyASqADne+56N+kzae/xx\n6NVLV/+ppk0bWzW0d2+tFyTpJ5ZO4JpYsxDefU3vcW1gc9T7NgN1ijm+xTue9r75xq4k+/VzHYmU\nx4MP2l7N48a5jkQkvuI137HAu8XNsKhB8qFQiFAoFM+vT6px46BzZ6hb13UkUh6VKtnmMRdeCJdc\nYiuIivhBOBwmHA6X+/NlGQXUAJgJnOE9zwVCwDaseWc+0IRIX8AI7342MBTY4L2nqXf8GuBS4OZi\nfittRgH9+CM0bGgdiaed5joaicXEiTB6tK3fVLmy62hEfi2Z8wBmAIWr2PcApkcdvxqoBDTEOntz\nsETxLdYfkAH8Ieozaev556FlSxX+6eCGG+z/44ABriMRiY/SZoqp2NX68Vh7/33Aa8A0oD7W2dsN\n2Om9fzDQC9gP9APe8o63BP4NVAZmAX1L+L20qAEUFNgSw6NGQdu2rqOReNi5E844w2oDbdq4jkbk\nl7QUhI+8956NHlm50saVS3qYM8f+vy5dqhnd4i9aCsJHnnwSbr1VhX+6advWZgn3Lan+KpIiVANI\nkK1bbRmB9es19j8d/fADtGgBw4fD//yP62hEjGoAPjF+vC0mpsI/PR19NDz7rNXwtm1zHY1I+agG\nkAD798NvfgOzZ1uHoaSve++FZcvgtdcgw6//miQwVAPwgTfftA1fVPinv6FDYeNGGxUkkmqUABJg\nwgQbJSLpr1IleO45GDgQ1q1zHY1I2fi10pqyTUDbt0OTJnZVWKWK62gkWUaOtCa/t9/WqC9xR01A\njj3/PHTposI/aO68E77/3jr/RVKFagBxVFBgC4X985+2raAEy/LlEArBJ59oL2FxQzUAh3JybA/Z\niy5yHYm4cNppcMcd8Oc/a+8ASQ1KAHE0YQL07KnhgEE2YIDtHvbvf7uOROTQ/FpUpVwT0N69ULs2\nLF6s6n/QffaZLRexZAnUCcSWR+IXagJyJDvbRv+o8JdmzeCWW+Dmm9UUJP6mBBAnU6fCNde4jkL8\n4t57YcMG20lMxK/UBBQHu3db88/q1XDiia6jEb/45BPo2NGahGrVch2NBIGagByYORNatVLhL7/U\nsiX06mULxon4kRJAHKj5R0oydKjND3j1VdeRiPyamoBitHOnrfy5caOWfpbivfsuXHutJQL9jUgi\nqQkoyWbOhNat9Q9bSnbJJdYXMGiQ60hEfikeCWAQsBxYBkwBjgBqAHOA1UA2UK3I+9cAuUC7OPy+\nU6++Cl27uo5C/G7kSNsz4P33XUciEhFrE1ADYB7QFPgJeBGYBZwGfA2MBAYC1YF7gCwsSZwD1AHe\nBk4B8ot8b0o0Ae3eDZmZtgxwjRquoxG/e+klGDbMJgsecYTraCQdJbsJ6FtgH3AUUNG7/wLoBEzy\n3jMJ6OI97gxM9T6zHsgDzo0xBmeys+Hss1X4S+n87nfQqJHVBkT8INYE8A3wKLARK/h3Yk0/NYHt\n3nu2e88BagOboz6/GasJpCQ1/0hZZGTAk0/CmDGQm+s6GhG7ao9FI+AOrCloF/AScH2R9xR4t5IU\n+9qwYcN+fhwKhQiFQuWPMgH27YPXX4eHHnIdiaSSevXgvvvgxhshHNbmMRKbcDhMOBwu9+dj7QPo\nDrQF/uQ9/wNwHnAZ0BrYBmQC84EmWD8AwAjvfjYwFFhU5Ht93wcwd66N6sjJcR2JpJoDB+CCC2zZ\n6D/96dDvFymtZPcB5GIFfmXvR9sAK4CZQA/vPT2A6d7jGcDVQCWgIdAYSMkiVM0/Ul4VKtjOYYMH\nw7ZtrqORIIvHRLABWCGfDyzGagNVgGlAfayztxvWPwAwGOgF7Af6AW8V852+rgEUFMBJJ9kcgNNP\ndx2NpKrBg2HtWnjxRdeRSLooaw1AM4HLYeVKuOIKW+1Rm79Iee3ZA2ecAU88YRPFRGKlmcBJMGuW\n/YNV4S+xqFwZxo6F226zOSUiyaYEUA5vvgkdOriOQtJBu3ZwzjkaTSZu+PUa1rdNQN99Z2v/b90K\nxxzjOhpJB198YbuIvfsuNG3qOhpJZWoCSrB582ztfxX+Ei+1a8Nf/wp9+mgLSUkuJYAyKmz/F4mn\nPn1g1y6YPNl1JBIkagIqg4ICW/v/rbdUVZf4y8mBzp1hxQqoXt11NJKK1ASUQLm5NvKnSRPXkUg6\nOvdcm1w4eLDrSCQolADKYO5caNNGwz8lcR5+GKZPh0VFF0cRSQAlgDKYOxcuu8x1FJLOqlWDRx6B\nm2+G/ftdRyPpTgmglA4cgHfeUQKQxLvuOusDePJJ15FIulMCKKUlS2z3r8xM15FIusvIsBnCDz5o\ncwREEkUJoJTmzdPVvyRPkyZw003wl7+4jkTSmRJAKc2dC5df7joKCZJ774WPP7ZhxyKJ4NfxLL6a\nB7B3Lxx3HGzcqPHZklyvvw533QVLl2ojeTk0zQNIgIULrUquwl+S7cor4ZRT4LHHXEci6UgJoBTe\neQd8tiWxBMiYMTBqFGza5DoSSTdKAKWwcCGcf77rKCSoTjoJbr3VmoJE4kl9AIdQUADHHw+ff64h\noOLOnj2QlWV7Cbdp4zoa8Sv1AcTZ6tVQpYoKf3GrcmV4/HG4/XYblCASD/FIANWA/wArgRVAK6AG\nMAdYDWR77yk0CFgD5ALt4vD7CfXhh2r+EX/o1AkaNrQ+AZF4iEcCGAPMApoCZ2IF+z1YAjgFmOs9\nB8gCunv37YGxcYohYRYuhPPOcx2FiM0QHjMG/vY32LLFdTSSDmItfKsCFwMTvOf7gV1AJ2CSd2wS\n0MV73BmYCuwD1gN5wLkxxpBQqgGInzRubDOE+/d3HYmkg1gTQEPgK2AisBgYDxwN1AS2e+/Z7j0H\nqA1sjvr8ZqBOjDEkzHffwdq10Ly560hEIgYPhvffh3DYdSSS6irG4fNnAbcBHwGPE2nuKVTg3UpS\n7GvDhg37+XEoFCLkYCB+To4V/pUqJf2nRUp09NEwejTcdpstUnj44a4jElfC4TDhGK4EYh0GWgv4\nEKsJAFyEdfKeBLQGtgGZwHygCZHkMMK7nw0MBYpuf+GLYaAPPgg7d9okHBE/KSiAK66ADh20YJxE\nJHsY6DZgE9bZC9AGWA7MBHp4x3oA073HM4CrgUpY0mgM5MQYQ8JoApj4VUYG/P3vtoPY1q2uo5FU\nFY+JYM2Ap7FCfS3QE6gATAPqY5293YCd3vsHA72wDuN+QHFrHTqvARQUwAkn2CJctWs7DUWkRPfc\nYyOCnnvOdSTiB2WtAWgmcAlWr4a2bWHDBqdhiBzU999D06YwZQpcfLHraMQ1zQSOEw3/lFRwzDHw\n6KM2Q/jAAdfRSKpRAiiBJoBJqvj9720z+fHjXUciqUYJoASqAUiqyMiAJ56AoUPhm29cRyOpRH0A\nxfjuO6hVC3bs0BwASR19+kCFCjY6SIJJfQBx8NFHmgAmqeeBB+DFF23pcpHSUAIohtr/JRUdd5w1\nA/XrZ8OYRQ5FCaAYav+XVHXTTfDVV/Dqq64jkVSgPoBf/TCceCJ8+inU8e0ydSIlmz8fevWCFSts\nIxkJDvUBxCgvz/7RqPCXVNW6NbRsqTWs5NCUAIr48EO1/0vqGzXKNo/ZtMl1JOJnSgBFaAE4SQcN\nGsCtt8KAAa4jET9TAijigw/gggtcRyESu4ED7e95wQLXkYhfqRM4yrffQmamJoBJ+njxRRg+HD75\nxCaJSXpTJ3AMcnLgrLNU+Ev66NYNqlaFp592HYn4kRJAFDX/SLqJXidoxw7X0YjfKAFEUQKQdNSs\nGXTtClHbbIsA6gP4WX4+1KhhG8GceGJSf1ok4b7+GrKyYN48OP1019FIoqgPoJxWrLAtIFX4Szo6\n/nj461/hjju0TpBEKAF41Pwj6e6WW2DbNpg+3XUk4hfxSgAVgCXATO95DWAOsBrIBqpFvXcQsAbI\nBdrF6fdjpgQg6a5iResQvusu+PFH19GIH8QrAfQDVgCFlct7sARwCjDXew6QBXT37tsDY+MYQ0yU\nACQILrvMOoXHjHEdifhBPArfukBH4GkinQ+dgEne40lAF+9xZ2AqsA9YD+QB58Yhhph8+aXdsrJc\nRyKSeI88Yrdt21xHIq7FIwE8BvQH8qOO1QS2e4+3e88BagObo963GXC+7ua778JFF2mmpATDySdD\nz54wZIjrSMS1ijF+/krgS6z9P1TCewqINA2V9PqvDIsatBwKhQiFSvr62IXDkMCvF/GdIUPg1FNh\nyRJo0cJ1NFJe4XCYcDhc7s/HOg/gYeAPwH7gSOBY4BXgHCwhbAMygflAEyJ9ASO8+9nAUGBRke9N\n6jyA00+HiRPhnHOS9pMizj31FEyebBdAGX6dESRlkux5AIOBekBD4GpgHpYQZgA9vPf0AAoHns3w\n3lfJ+0xjICfGGGLy1Ve2ZrqugiRoeveGnTvhlVdcRyKuxHsETuFl+wigLTYM9DIiV/wrgGne/ZtA\nHw7ePJRwhe3/FWNtDBNJMRUqwGOPQf/+GhYaVH6t+CWtCej226FePW2cIcHVtSu0agX33HPo94q/\nlbUJKPAJ4Mwzbancc50PRhVxIy/PtkH9/HOoVct1NBILJYAy+PpraNQI/u//1AQkwda/vy0XrX0D\nUpsWgyuDd9+12b8q/CXohgyBN96wYaESHIFOAHPmQJs2rqMQca9qVfjf/9VqoUET6ASQnQ3tfLMc\nnYhbGhYaPIFNAGvXwu7d2hxDpJCGhQZPYBNA4dW/ZkCKRBSuFvr4464jkWQIdAK44grXUYj4zyOP\nwKhRsHWr60gk0fx6/ZvQYaD79tn2j9r/V6R4/fvDN9/AM8+4jkTKQsNAS2HRIjjpJBX+IiUpHBa6\neLHrSCSRApkAsrOhbVvXUYj4V9WqcP/9Ghaa7gKZAF5/Ha680nUUIv7Wuzfs2gUvv+w6EkmUwPUB\nFC79vG2bZgCLHMr8+dCrF6xcCUce6ToaORT1ARzCzJnQsaMKf5HSaN0amje3+QGSfgKXAGbMgE6d\nXEchkjo0LDR9BaoJ6LvvoE4d2LIFqlSJ+9eLpK2777b+gPHjXUciB6MmoIPIzobzz1fhL1JWQ4ZY\n7fmzz1xHIvEUqATw2mtw1VWuoxBJPdWqwX33wV13aVhoOglME9BPP0Fmpu16VLt2XL9aJBD27bMd\n9EaNgt/+1nU0UpxkNwHVA+YDy4HPgb7e8RrAHGxT+GygWtRnBgFrgFwgaYsxZ2fDGWeo8Bcpr8MP\nt8L/7rstGUjqizUB7AP+ApwGnAfcCjQF7sESwCnAXO85QBbQ3btvD4yNQwyl8uKL0K1bMn5JJH11\n7Ah168JTT7mOROIh3k1A04F/eLdLge1ALSAMNMGu/vOBv3nvnw0MAxYW+Z64NgHt2WPNP7m52vRa\nJFZLl9pSKqtWWd+A+IfLUUANgBbAIqAmVvjj3df0HtcGNkd9ZjNQJ44xFGv2bDjrLBX+IvFw5pk2\nl+ahh1xHIrGK13zYY4CXgX7Ad0VeK/BuJSn2tWHDhv38OBQKEQqFyh3cCy9A9+7l/riIFPHAA7ab\n3s03Q6NGrqMJrnA4TDgcLvfn49EEdDjwOvAmULiPUC4QArYBmVhHcRMifQEjvPvZwFCs1hAtbk1A\nO3ZAw4bw3/9CjRpx+UoRwWoAn34KL73kOhIplOwmoAzgGWAFkcIfYAbQw3vcA+sbKDx+NVAJaAg0\nBnJijOGgXnjBdv5S4S8SX3feaXtrvPee60ikvGJNABcC1wOtgSXerT12hd8WGwZ6GZEr/hXANO/+\nTaAPB28eitnEidCzZyJ/QSSYKleG4cMtEeTnu45GyiOtJ4J9/jm0bw8bNkCFCnGISkR+IT8fzjsP\n+vWD665zHY1oLaAoEyfCH/+owl8kUQ47DEaPhkGDYPdu19FIWaVtDWD3bqhfH3JybP9fEUmcbt1s\neOiQIa4jCbay1gDSNgGMH2+bv8yYEaeIRKRE//0vnHOONbtmZrqOJriUALDVCps1g0cf1ebvIsky\nYAB88w08/bTrSIJLfQDAO+/A/v3Qpo3rSESCY/Bgq3Vrz4DUkZYJYNQo6NsXMvxavxFJQ9WqwdCh\nNixUewakhrRLAEuW2O2GG1xHIhI8N95oewe/8YbrSKQ00i4BPPSQrVd+5JGuIxEJnooVtWdAKkmr\nBLB8OSxYYFchIuJGhw42BPtf/3IdiRyKX1vJyzUKqFMnuPRS27dURNxZtswGYeTmQvXqrqMJjsAO\nA50/H3r3hpUr4YgjEhSViJTajTfCscdak5AkRyATQH6+TUIZMEDr/ov4xbZttmfAokXaMyBZAjkP\nYNw4W5lQe/6K+EetWjYkdOBA15FISVK+BrBpE7RoYZ2/TZsmOCoRKZM9e6BJE3j+ebj4YtfRpL9A\n1QDy862dsW9fFf4iflS5MowYoT0D/CqlE8CoUbBrly1FKyL+dPXVtiT7lCmuI5GiUrYJ6J13rM3/\no49szLGI+NcHH9gAjVWr4KijXEeTvgLRBLR8uRX+kyer8BdJBRdcABdeCCNHuo5EorlKAO2BXGAN\nUKYxArm5ts3j6NFa7VMklYwcCf/4B2zc6DoSKeQiAVQA/oElgSzgGqBUXbgffACtW9t6P9p/VCS1\n1K8Pt98O/fu7jkQKuUgA5wJ5wHpgH/AC0PlgH9i7F4YPh65d4ZlnbJ9fEUk9/fvDwoXWhyfuuUgA\ndYBNUc83e8d+5fvvrcA//XR4913r8O3YMSkxikgCHHUUPPII9OsHBw64jkYqOvjNUs3wOv/8yIJS\nY8eqvV8kXfz+9/Dkk1ajr1XLdTTB5iIBbAHqRT2vh9UCfuHUU4fRujVUqgQVK4aAUHKiE5GEysiw\nOQHaNCZ2q1aFWb06XO7Pu5gHUBFYBVwOfAHkYB3BK6PeE9Om8CIiQVTWeQAuagD7gduAt7ARQc/w\ny8JfRESSIGVnAouIyC8FYiawiIjETglARCSglABERAJKCUBEJKCUAEREAkoJQEQkoJQAREQCSglA\nRCSglABERAJKCUBEJKCUAEREAkoJQEQkoJQAREQCSglARCSglABERAJKCUBEJKCUAEREAkoJQEQk\noJQAREQCKpYE8Ai2mftnwCtA1ajXBgFrgFygXdTxlsAy77UxMfy2iIjEKJYEkA2cBjQDVmOFPkAW\n0N27bw+MJbJJ8TigN9DYu7WP4fcDIRwOuw7BN3QuInQuInQuyi+WBDAHyPceLwLqeo87A1OBfcB6\nIA9oBWQCVYAc733PAl1i+P1A0B93hM5FhM5FhM5F+cWrD6AXMMt7XBvYHPXaZqBOMce3eMdFRMSB\niod4fQ5Qq5jjg4GZ3uN7gb3AlDjGJSIiPncD8D5wZNSxe7xbodlYE1AtrNO40DXAP0v43jygQDfd\ndNNNtzLd8kiS9sBy4Pgix7OAT4FKQENgLZFO4EVYMsjAmozUCSwikoLWABuAJd5tbNRrg7FMlAtc\nEXW8cBhoHvBEcsIUERERERHfao/VGtYAAx3HkmwTgO1YDalQDawjfjU276Kag7hcqAfMx5oYPwf6\neseDeD6OxJpOPwVWAMO940E8F4UqYK0OhQNRgnou1gNLsXNROLw+Zc9FBaxpqAFwOPYH39RlQEl2\nMdCCXyZu1KSUAAAB/0lEQVSAkcAA7/FAYESyg3KkFtDce3wMsAr7Wwjq+TjKu68ILAQuIrjnAuBO\nYDIww3se1HOxDivwo6XsuTgfGzFUqOhooiBowC8TQC5Q03tcy3seRNOBNuh8HAV8hM3AD+q5qAu8\nDbQmUgMI6rlYBxxX5FiZzoWfFoOrA2yKel44gSzIamLNQnj3NQ/y3nTVAKsZLSK45+MwrEa8nUjT\nWFDPxWNAfyKrEEBwz0UBlgw/Bv7sHSvTuTjURLBkKnAdgM8VjvMNkmOAl4F+wHdFXgvS+cjHmsSq\nAm9hV7/RgnIurgS+xNq8QyW8JyjnAuBCYCtwAtbuX/Rq/5Dnwk81gC1Y51+hevxy6Ygg2k5kJnYm\n9scfFIdjhf9zWBMQBPt8AOwC3sCGUwfxXFwAdMKaPqYCl2F/H0E8F2CFP8BXwKvAuZTxXPgpAXyM\nrRDaAJtE1p1IJ09QzQB6eI97ECkI010G8Aw26uXxqONBPB/HExnJURloi10BB/FcDMYuDBsCVwPz\ngD8QzHNxFLa4JsDR2LL7y0jxc9EBG/GRR2R56aCYCnyBrau0CeiJ9fC/TQoO6YrRRVizx6dEJhq2\nJ5jn4wxgMXYulmLt3xDMcxHtUiIXiEE8Fw2xv4lPsaHSheVlEM+FiIiIiIiIiIiIiIiIiIiIiIiI\niIiIiIiIiEj6+n8kJLxlzmZuSwAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0xb08891ec>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEACAYAAAC6d6FnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHVdJREFUeJzt3XmUFNXZx/HvMICyKSI6LALDYZ+EsAmooLYoiAYFNAom\nIkYkMSaYmMTXjbxMyEk0JjHKMZiTBBQNwSUmwEDYHGghMSyDrI7DFkGGXVEDr6gD9PvHrbabZoCZ\n6eq+tfw+59SpmuqlHutgPV33uXUviIiIiIiIiIiIiIiIiIiIiIiIiIhICLUClgBvAxuB+5z9hUA5\nsMZZrkv6zMPAFqAMGJStQEVExF3NgO7OdkNgE9AFmAD8sJL3FwBrgTpAPrAVqJXxKEVE5CTpXnz3\nYi7oAIeBd4CWzt85lbx/KDADqAC2YxJAnzRjEBGRGnDz13c+0ANY7vw9DlgHTAEaO/taYJqG4spJ\nJAwREckitxJAQ+CvwPcxdwLPAm0xzUN7gN+c5rMxl2IQEZFqqO3Cd9QBXgP+DMx09u1Pev1PQJGz\nvQtTOI67yNl3gnbt2sW2bdvmQmgiIqGyDWhf1TeneweQg2niKQWeStrfPGl7OLDB2Z4NjATqYu4Q\nOgArU79027ZtxGIxLbEYEyZMsB6DVxadC50LnYvTL0C76lzA070D6AfcDqzHdPcEeAS4DdP8EwPe\nBb7tvFYKvOKsjwL3oiYgEREr0k0A/6Tyu4h5p/nML5xFREQsUh98j4tEIrZD8AydiwSdiwSdi5qr\nrK++F8Sc9iwREaminJwcqMZ1XXcAIiIhpQQgIhJSSgAiIiGlBCAiElJKAJUoL4c77rAdhYhIZikB\nVKKoCF58EdQRSUSCTAmgEjt2mPXBg3bjEBHJJCWASqxebdZbttiNQ0Qkk5QAUsRiJgFEIrB1q+1o\nREQyRwkgxfbtUK8e9O+vOwARCTYlgBSrV0OvXtChg+4ARCTYlABSxBNA+/a6AxCRYFMCSKE7ABEJ\nCyWAJPECcK9e0LQpHDsGH3xgOyoRkcxQAkiyYwecdRY0bw45OboLEJFgUwJIEv/1H9e+vRKAiASX\nEkCS1ATQoQNs3mwvHhGRTFICSJKaALp0gbIye/GIiGSSEoAjuQAcV1AApaX2YhIRySQlAMd770Gd\nOtCiRWJfp06mBnD0qL24REQyRQnAkfrrH8yQEC1bwrZtdmISEcmkdBNAK2AJ8DawEbjP2d8EWARs\nBhYCjZM+8zCwBSgDBqV5fNdUlgBAzUAiElzpJoAK4H7gS8AlwHeBLsBDmATQESh2/gYoAEY468HA\nZBdicIUSgIiETboX373AWmf7MPAO0BK4EZjm7J8GDHO2hwIzMIljO7AV6JNmDGmrrAAcpwQgIkHl\n5q/vfKAHsALIA/Y5+/c5fwO0AMqTPlOOSRhW7dwJubknFoDjlABEJKhqu/Q9DYHXgO8Dh1JeiznL\nqVT6WmFh4RfbkUiESCSSVoCnE//1n5Nz8mudO8OmTWZcoNzcjIUgIlJt0WiUaDRa489XcsmrtjrA\nHGAe8JSzrwyIYJqImmMKxZ1J1AIed9bzgQmYu4ZksVgWZ2QfPx5q1YKJEyt/PT8fiouhXbushSQi\nUm055ldsla/r6TYB5QBTgFISF3+A2cBoZ3s0MDNp/0igLtAW6ACsTDOGtL31FvTseerXv/xlWL8+\ne/GIiGRDugmgH3A7cBWwxlkGY37hD8R0Ax1A4hd/KfCKs54H3Mvpm4cy7nQF4Lhu3WDduuzFJCKS\nDW40AWVC1pqAysvNr/99+yqvAQC8+ipMnw4zZ1b+uoiIF2S7Ccj3TlcAjuveXXcAIhI8SgBnaP4B\nU/x9/3346KPsxCQikg1KAFVIALVqwVe+orsAEQmWUCeAqhSA47p1g7Vrz/w+ERG/CHUC2L3bPODV\nqtWZ39u9uxKAiARLqBNAVQrAcUoAIhI0SgBVaP4B8zDYpk3w+eeZjUlEJFuUAKqYAOrXh/bt9USw\niASHEkAVEwBA376wInXUIhERnwptAti9GyoqoHXrqn+mTx9YaX3kIhERd4Q2AVSnABynOwARCZLQ\nJ4DqKCgwYwfpiWARCQIlgGqoXdsMHFdSkpmYRESySQmgmtQMJCJBEcoEsGcPfPYZtGlT/c+qECwi\nQRHKBFCTAnBc377w73+bcYRERPws1AmgJlq3hnr1YPNmd2MSEck2JYAauPxyWLbMvXhERGxQAqgB\nJQARCYLQJYC9e+HIEcjPr/l3XHEFLF3qWkgiIlaELgGsXm368tekABzXuTMcPgw7d7oXl4hItoUy\nAaTT/AMmeagZSET8TgmghpQARMTv3EgAU4F9wIakfYVAObDGWa5Leu1hYAtQBgxy4fjVsno19O6d\n/vdccQVEo+l/j4iILW4kgOeAwSn7YsCTQA9nmefsLwBGOOvBwGSXYqiSPXvSLwDHde8OBw6YweFE\nRPzIjYvvMuDDSvZXVmYdCswAKoDtwFagjwsxVMnq1XDxxekVgONyc2HAAHj99fS/S0TEhkz++h4H\nrAOmAI2dfS0wTUNx5UDLDMZwgpISkwDcMmgQLFzo3veJiGRT7Qx977PARGf7Z8BvgDGneG+lo+oU\nFhZ+sR2JRIhEImkHVVICY04VRQ0MHAiPPALHj0Ot0JXTRcS2aDRKNI1ipAuNIQDkA0VA1zO89pCz\n73FnPR+YAKQOsByLuTzaWiwGzZvDqlXQqpV739uxI7zyiqkJiIjYlGPat6t8Xc/U79bmSdvDSfQQ\nmg2MBOoCbYEOQFYGV961yySBiy5y93vVDCQifuVGApgBvAl0AnYCdwG/BNZjagBXAvc77y0FXnHW\n84B7OUUTkNtKSmo+BPTpDByoBCAi/uTy5dA1rjcBjR9v2uknTjzze6vj0CFo2dLcYTRq5O53i4hU\nh1eagDzH7R5AcY0awWWX6S5ARPwnFAkgFstcAgC44QaYPTsz3y0ikimhSAA7dkDdutCiRWa+/4Yb\nYO5cOHo0M98vIpIJoUgAmfz1D2aayFatzFzBIiJ+oQTgkhtugKKizB5DRMRNSgAuufFGmDXL1BtE\nRPwg8AkgFnNvDoDT6dULPv0U3n47s8cREXFL4BPAf/5jumrm5WX2ODk5cMst8PLLmT2OiIhbAp8A\nstH8EzdihEkAagYSET9QAnDRxRebrqDr1mXneCIi6VACcFFODtx6q5qBRMQfAj0W0PHjcN55pg5w\n/vkuRFUFa9bAzTfDtm3uDzwnInI6GgsoyebN0KRJ9i7+YOYFqF0bVqTOcCAi4jGBTgArVkDfvtk9\nZk4OjB4N06Zl97giItWlBJABd9xhZgk7ciT7xxYRqapAJ4CVK6FPn+wft1UrU3ieNSv7xxYRqarA\nJoAjR6C0FHr2tHP8O++E556zc2wRkaoIbAJYswa6dIF69ewcf9gw0wW1vNzO8UVEziSwCWDlSjvt\n/3H16plnAnQXICJeFdgEsGKFnfb/ZPfcA3/4gyaKERFvCnQCsHkHANCtG7Rpo+kiRcSbApkADhyA\ngwehUyfbkcB3vwu/+53tKEREThbIBLBypemGWcsD/3U332x6I73zju1IRERO5MYlciqwD9iQtK8J\nsAjYDCwEGie99jCwBSgDBrlw/JN4ofknrm5duPtumDzZdiQiIidyIwE8BwxO2fcQJgF0BIqdvwEK\ngBHOejAw2aUYTuClBADw7W/D9Onw0Ue2IxERSXDj4rsM+DBl341AfDScacAwZ3soMAOoALYDWwFX\n++rEYrBqlf0eQMkuugiGDIFnn7UdiYhIQqZayfMwzUI46/iEjC2A5EejyoGWbh54yxYzBWSzZm5+\na/oeeAAmTTLzBouIeEHtLBwj5iyne/0khYWFX2xHIhEikUiVDua15p+4rl2hRw948UUYO9Z2NCIS\nBNFolGg0WuPPuzVlST5QBHR1/i4DIsBeoDmwBOhMohbwuLOeD0wAUkfPr/GEMOPGQX4+/OhHNfp4\nRkWjph5QWgq5ubajEZGg8cqEMLOB0c72aGBm0v6RQF2gLdABWOnmgb3wBPCpXHklnHsuzJx55veK\niGSaG3cAM4ArgaaY9v7/BWYBrwCtMcXeW4F4H5hHgLuAo8D3gQWVfGeN7gCOHIGmTc2DYPXrV/vj\nWVFUBOPHm8HqvPCcgogER3XvALw6a22NEsCyZabpZ6Wr9xTuisXMHcqDD8LXvmY7GhEJEq80AVnx\n5pvQr5/tKE4vJwcmToQJE+DYMdvRiEiYBSoB/OtfcNlltqM4s8GD4Zxz4NVXbUciImEWmCagWAwu\nuADWrjUPXnndokWmx9LGjVA7G51xRSTwQtsEtGULNGjgj4s/wDXXQPPmMHWq7UhEJKwCkwDefNMf\nzT9xOTnw61+bWsChQ7ajEZEwClQC8HoBOFWvXuZO4IknbEciImEUmBrAl78ML7wAPXtmKKIM2bkT\nuneHdev803wlIt4UyucAPvwQWrc2az8WVB99FHbtgueftx2JiPhZKIvAy5dD797+vPiDeShs0SL4\n979tRyIiYRKIBOC3AnCqc84xBeHvfAeOHrUdjYiERWASgN8KwKlGjoTzz9cE8iKSPb6vARw9Ck2a\nwI4dcN55GY4qw8rKoH9/WL8eWrSwHY2I+E3oagDr10OrVv6/+AN07mzmC/jhD21HIiJh4PsE8M9/\n+r/5J9mjj0JJCcyaZTsSEQk63yeApUvNRCtBUb++GR7i3nvh4EHb0YhIkPm6BhCLQV4erF5tmoGC\n5L774KOPzMNtIiJVEaoaQFmZGQAuaBd/gMceM8Nbz5ljOxIRCSpfJ4ClS+GKK2xHkRkNGsCUKXDP\nPfDBB7ajEZEg8n0CCFL7f6pIBG69Fb71LdPcJSLiJt8mgFgM3ngjuHcAcY89Btu2wR//aDsSEQka\n3xaB333XdP/ctcuMrR9k77wDl19uJr3v0sV2NCLiVaEpAsfb/4N+8Qdz0f/FL+C22+Czz2xHIyJB\n4fsEEBZjx0K7dvDjH9uORESCItMJYDuwHlgDrHT2NQEWAZuBhUDjmnxx2BJATo7pFTRvHkyfbjsa\nEQmCTDegvAv0ApKfaX0CeN9ZPwicBzyU8rnT1gB274auXeHAAajl23uYmlm3zkwjuWSJmQVNRCTO\nizWA1GBuBKY529OAYdX9wqVLTVE0bBd/gG7d4Mkn4aab4OOPbUcjIn6W6UtoDHgdKAHGOvvygH3O\n9j7n72qJRsPV/JNq1ChzF3DnnXo+QERqLtOTKPYD9gAXYNr9y1JejznLSQoLC7/YjkQiRCKRL/5e\nvNjMnhVmv/2teVBs4kSYMMF2NCJiQzQaJRqN1vjz2exEOQE4jLkTiAB7gebAEqBzyntPWQPYuRN6\n9oR9+8LZBJRs717o2xeeeAJGjLAdjYjY5qUaQH2gkbPdABgEbABmA6Od/aOBmdX50uJiuOoqXfwB\nmjWD2bNh3DhYufLM7xcRSZbJJqA84O9Jx5mO6fZZArwCjMF0E721Ol+6eDFcfbV7Qfpdt27wpz/B\n8OGwYgVcdJHtiETEL7z6HG2lTUCxmLnARaPQoUP2g/KyJ56AGTPMcBENG9qORkRs8FITkOs2bYLc\nXGjf3nYk3vPAA9CrF9xyC1RU2I5GRPzAVwkg3vwThvF/qisnB37/e5Mgx45V91AROTNfJYDiYhgw\nwHYU3lW7Nrz8spkp7dFHbUcjIl7n1d/SJ9UAjh2DCy+E9euhZUtLUfnE+++bobLHjYPvfc92NCKS\nLdWtAWT6QTDXrFsHF1ygi39VNG0K8+dD//6Ql2fqAiIiqXyTAIqL1f2zOtq2hblz4dprzfzC119v\nOyIR8Rrf1ABef10JoLq6d4dZs2D0aDN6qIhIMl/UAD75xDRllJfDuedajMqnolHTDDR7Nlx6qe1o\nRCRTAvkcwNKl5tesLv41E4nAtGkwbBisWWM7GhHxCl8kgPnzYfBg21H42/XXw+TJZl1aajsaEfEC\nXxSB58/XNIhuuPlm05w2cKCpqXTpYjsiEbHJ8wlg+3Y4eBB69LAdSTCMGmWeEr76apMECgpsRyQi\ntng+ASxYAIMGafhnN91xhzmf11wDixbBl75kOyIRscEXCeCmm2xHETy3326SwMCBsHChJpgXCSNP\ndwOtqDBP/27ebIaBEPe99BLcf79JAl272o5GRNIRqKEgli+Hdu108c+kkSPNSKKDBplie7dutiMS\nkWzxdAJQ98/sGDHCDCN97bXmYbE+fWxHJCLZ4OnS6ty5SgDZ8rWvmaklhwyBN96wHY2IZINnawA7\ndsTo2RP27jXj3Et2LF5smoVeeEHJV8RvAjMUxJw55qlVXfyza8CAxAByr71mOxoRySTPXl6LiuCu\nu2xHEU6XXmrqL9dfb54cHjXKdkQikgmebQJq1ChGeTmcc47tUMKrtNT0Dho/Hu65x3Y0InImgekG\n2revLv62FRSYgvA118ChQ/DAA7YjEhE32aoBDAbKgC3Ag5W94cYbsxqPnEK7drBsGUyZAj/5iRlH\nSESCwUYTUC6wCbgG2AWsAm4D3kl6T+zdd2Pk52c/OKnc/v3mOYF+/WDSJI3NJOJFfugF1AfYCmwH\nKoCXgKGpb9LF31suvNDMLLZ+vSkKV1TYjkhE0mUjAbQEdib9Xe7sE48791wzON9//2tmF/vkE9sR\niUg6bBSBq9SKXFhY+MV2JBIhEolkKBypjnr14G9/M110r73WdNdt3Nh2VCLhFI1GiUajNf68jRrA\nJUAhphAM8DBwHPhl0ntOmBRevOf4cfjBD8x8zQsWQF6e7YhExA81gBKgA5AP1AVGALMtxCFpqFUL\nnn4ahg+H/v3NzG0i4i82moCOAt8DFmB6BE3hxB5A4hM5OTBhAjRpApdfbu4ENMWkiH949klgNQH5\ny/Tp8KMfmXGE+va1HY1IOFW3CUgJQFwzZw5885swY4Z5elhEsssPNQAJqCFDzAiiX/+6RhIV8QPP\njgUk/nTFFaYW8NWvwkcfwZgxtiMSkVNRAhDX9ehhBpEbNAgOHtQgciJepRqAZEx5uXlYbMgQePxx\n02tIRDJHRWDxlA8+MM1BHTuaOYfr1rUdkUhwqQgsnnL++Wae4f/+18ww9vHHtiMSkTglAMm4+vVN\nr6BOnUyReNcu2xGJCCgBSJbk5sIzz8A3vgGXXQYbN9qOSERUA5Cs+8tf4P774aWX4KqrbEcjEhyq\nAYjnff3r5uI/cqRJBiJih+4AxJqNG00X0VGj4Kc/1TSTIulSN1Dxlf374aaboFkzmDYNGjSwHZGI\nf6kJSHzlwguhuBgaNTLzCuzceebPiIg7lADEurPOgqlT4fbb4ZJLYPly2xGJhIOagMRT5s41Q0r/\n6lcwerTtaET8RTUA8b3SUjPV5IAB8NRT5g5BRM5MNQDxvYICWLUKDhwwU02+957tiESCSQlAPOmc\nc+DVV2HECOjTBxYutB2RSPCoCUg87403zMNj3/oWjB9vhpUQkZOpBiCBtHu3eWCsogL+/Gdo3dp2\nRCLeoxqABFKLFrBokZlb4OKL4a9/tR2RiP/pDkB8Z9Uq0yR05ZXw9NN6elgkzit3AIVAObDGWa5L\neu1hYAtQBgzK0PElwHr3hrfegqNHoVs3WLrUdkQi/pSpO4AJwCHgyZT9BcBfgN5AS+B1oCNwPOV9\nugOQKpk1C+69F26+GR57THcDEm5euQM4VRBDgRlABbAd2Ar0yWAMEnBDh8KGDWaqya5dYckS2xGJ\n+EcmE8A4YB0wBWjs7GuBaRqKK8fcCYjUWJMmZiTRSZNMT6GxY81k9CJyerXT+OwioFkl+x8FngUm\nOn//DPgNMOYU31NpW09hYeEX25FIhEgkUsMwJSyGDIG334af/MQ8Tfzzn8Ndd2meAQmuaDRKNBqt\n8eez0QsoHygCugIPOfsed9bzMfWCFSmfUQ1A0rJmjakNAEyeDD162I1HJBu8UgNonrQ9HNjgbM8G\nRgJ1gbZAB2BlhmKQEOvRA/71L7j7brjuOhgzBnbtsh2ViLdkKgH8EliPqQFcCdzv7C8FXnHW84B7\nOUUTkEi6atUyF/6yMrjgAvjKV0zz0KFDtiMT8QY9CCah8d57ZiyhRYvM+u67NdS0BItXmoBEPKd1\na3jhBTPpzLx50L49PPMMfPqp7chE7FACkNDp2RPmzIG//Q0WLDCJYNIkOHLEdmQi2aUEIKHVuzcU\nFZmniRcvhjZtTI1g717bkYlkhxKAhF6vXjBzJixbZh4g69IF7rwT1q2zHZlIZqkILJLi4EH4wx9M\nfaB1a/Nk8a23apwh8T5NCCPikqNH4R//gD/+0TxTMGKE6TnUsyfkePX/HAk1JQCRDCgvh+efh6lT\noW5duO02GDkSOnWyHZlIghKASAbFYrBiBbz0Erz8MjRvbu4Mhg1TMhD7lABEsuTYMTMZzcsvm95E\nDRrADTeYpV8/qFPHdoQSNkoAIhbEYmYAuqIis/znPzBgQGLp1El1A8k8JQARD9i9G4qLzfMFxcXm\nbmHAALjqKrj0UpMQNEy1uE0JQMRjYjFzR1BcDG+8AcuXm66mffrAJZeYpXdvaNrUdqTid0oAIj6w\nf78pJi9fbpaSEmjUyIxY2q1bYunQAWqnM22ThIoSgIgPxWKwfbt5+jh52bMH2rUzTUYdOybWHTvC\n+efbjlq8RglAJEAOH4YtW2Dz5sSyaZNZatc24xe1aWOeWE5et2lj5kBQ4TlclABEQiAWg/ffhx07\nzDwHO3acvH3oEFx4ITRrZpa8vMR2/O8mTRKL5kbwPyUAEQHMPAf795vRTffuhX37Etvxvw8eTCx1\n6sB555lkEF83aQKNG5v6RKNG0LDhievUfQ0bQm6u7f/y8KpuAlB5SSSgzj7bNAm1bn3m98Zi8Mkn\nJyaEDz9MrA8fNnMqHz5s7iwOHUpsJ+/7v/8zieTss0+/nHWWWdepo2Yqm7x66nUHIOJDx4/D55+b\nu49PP4XPPktsV7ZUVNiOOFhGjVITkIhIKGlOYBERqRIlABGRkEonAdwCvA0cA3qmvPYwsAUoAwYl\n7e8FbHBeezqNY4uISJrSSQAbgOHA0pT9BcAIZz0YmEyiTepZYAzQwVkGp3H8UIhGo7ZD8AydiwSd\niwSdi5pLJwGUAZsr2T8UmAFUANuBrUBfoDnQCFjpvO8FYFgaxw8F/eNO0LlI0LlI0LmouUzUAFoA\n5Ul/lwMtK9m/y9kvIiIWnOlBsEVAs0r2PwIUuR+OiIj4yRJOLAI/5Cxx8zFNQM2Ad5L23wb8/hTf\nuRWIadGiRYuWai1bybIlmN49cQXAWqAu0BbYRqIIvAKTDHKAf6AisIiILw0HdgJHgL3AvKTXHsFk\nojLg2qT98W6gW4FJ2QlTREREREQ8azDmrmEL8KDlWLJtKrAPc4cU1wRTiN8MLAQaW4jLhlaYpsW3\ngY3Afc7+MJ6PszFNp2uBUuAxZ38Yz0VcLrCGREeUsJ6L7cB6zLmId6/37bnIxTQN5QN1MP/gu9gM\nKMsuB3pwYgJ4AvgfZ/tB4PFsB2VJM6C7s90Q2IT5txDW81HfWdcGlgP9Ce+5APghMB2Y7fwd1nPx\nLuaCn8y35+JSTI+huNTeRGGQz4kJoAzIc7abOX+H0UzgGnQ+6gOrgC8R3nNxEfA6cBWJO4Cwnot3\ngdSZoat1Lrw0GFxLTFE5Lv4AWZjlYZqFcNZ5p3lvUOVj7oxWEN7zUQtzR7yPRNNYWM/Fb4EHgONJ\n+8J6LmKYZFgCjHX2VetceGlGsJjtADwu3s83TBoCrwHfBw6lvBam83Ec0yR2LrAA8+s3WVjOxRBg\nP6bNO3KK94TlXAD0A/YAF2Da/VN/7Z/xXHjpDmAXpvgX14oTh44Io30knsRujvnHHxZ1MBf/FzFN\nQBDu8wHwMTAX0506jOfiMuBGTNPHDGAA5t9HGM8FmIs/wAHg70AfqnkuvJQASjAjhOZjHiIbQaLI\nE1azgdHO9mgSF8KgywGmYHq9PJW0P4znoymJnhz1gIGYX8BhPBePYH4YtgVGAouBUYTzXNTHDK4J\n0AAz7P4GfH4ursP0+NiKmVMgTGYAu4HPMbWQb2Iq/K/jwy5daeqPafZYi7nYrcF0EQ7j+egKvIU5\nF+sx7d8QznOR7EoSPxDDeC7aYv5NrMV0lY5fL8N4LkRERERERERERERERERERERERERERERERERE\nguv/AVirL6unWV+kAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0xb050372c>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEACAYAAAC9Gb03AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEXNJREFUeJzt3W2MXNV9x/HvetdPa9debY3WC956XTuATYwLJYE0RJkQ\nSN2KEKJIxZGIXCrlTRtBqZRg4AX7KqWpqlAlyos2pHWqhCZqXGTaPNgB5kWV8qRgMMYbYwswLHj9\ngB9SBMY22xfnjmd28a53987sfTjfjzSae+/OzD0cVj8f/8+ZY5AkSZIkSZIkSZIkSZIkSdI0fA8Y\nBnY2XOsGtgN7gG1AV8PP7gZeAgaBz8xQGyVJKXwCuILRQf8N4GvJ8V3A/cnxGmAHMBvoB/YCs2ak\nlZKkVPoZHfSDQE9yvDQ5hzCav6vhdT8Hrml14yRJ55ZmpN1DKOeQPNdC/0Lg9YbXvQ5clOI+kqQU\nmlVSGUkeE/1ckpSBjhTvHSaUbA4AvcDB5PoQ0NfwumXJtVHa21eOnDmzL8XtJSlK+4BVU3lDmhH9\nVmBjcrwReLjh+gZgDrAC+BDw1Ng3nzmzjy1bRhgZ8XHfffdl3oa8POwL+8K+mPgBrJxqWE826B8C\nfgVcArwG3EZYZXMDYXnlddRX3bwI/Dh5/hnwl4xTunn++ak2V5I0VZMt3XxxnOvXj3P968ljQga9\nJLVepuvbn3suy7vnR6VSyboJuWFf1NkXdfZFOm0Z3nuks3OE4WFYuDDDVkhSgbS1tcEUszvTEf3q\n1fDCC1m2QJLKL9Ogv/xy6/SS1GqZB711eklqrUyDft06R/SS1GqZTsYeOjTCypVw7Bi0ZdkSSSqI\nwk3GLlkSVtzs359lKySp3DLfJ94JWUlqrVwEvROyktQ6mQe9E7KS1FqZB72lG0lqrUxX3YyMjHDq\nFCxaBEeOQGdnhq2RpAIo3KobgNmz4ZJLYNeurFsiSeWUedCDdXpJaqVcBL11eklqHYNekkou88lY\ngOFhWLMGDh92KwRJmkghJ2MBenqgowOGhrJuiSSVTy6CHpyQlaRWyU3QW6eXpNYw6CWp5HIV9G5u\nJknNl4tVNwAnT0JXFxw9CvPmZdgqScqxwq66AZg7F1atgt27s26JJJVLboIerNNLUisY9JJUcrkL\neidkJam5mhH0dwO7gJ3AD4G5QDewHdgDbAO6JvNB69aFoG+Yo5UkpZQ26PuBLwNXAmuBdmADsIkQ\n9BcDjybn59XbG0J+eDhlqyRJZ6UN+hPAKaAT6Eie3wBuAjYnr9kM3DyZD2trs04vSc2WNujfAv4B\n2E8I+GOEkXwPUBuXDyfnk2KdXpKaK23QrwT+mlDCuRBYCNw65jUjyWNS3NxMkpqrI+X7rwJ+BRxJ\nzrcAHwMOAEuT517g4LnePDAwcPa4UqlQqVS4/HJ44IGUrZKkkqhWq1Sr1VSfkXYLhHXAD4CPAO8C\n/wo8BSwnhP/fESZiu/jghOyoLRBq3nkHurvh+HGYMydl6ySpZLLYAuE54PvAM0Ct4PJPwP3ADYTl\nldcl55Myfz7098PgYMqWSZKAHG1q1uiWW+Czn4Vbx1b7JSlyhd7UrJETspLUPLkMetfSS1LzGPSS\nVHK5DPq+vrD65tChrFsiScWXy6B3KwRJap5cBj0Y9JLULAa9JJVcroPezc0kKb1cfmEK4O234YIL\n4MQJ6Ei7I48klURpvjAFsGABLFsGe/Zk3RJJKrbcBj1Yp5ekZjDoJankch/0TshKUjq5Dno3N5Ok\n9HId9MuXh3+A5K23sm6JJBVXroN+1ixYuxZ27sy6JZJUXLkOerBOL0lpFSLordNL0vTlPuidkJWk\ndHK7BULNiRPQ2xue29tnoFWSlGOl2gKhZtEi6OmBffuybokkFVPugx6ckJWkNAoR9NbpJWn6ChH0\nrryRpOkz6CWp5HK/6gbgzBlYvBiGhsKzJMWqlKtuICyrvOwyt0KQpOkoRNCDE7KSNF2FCXrr9JI0\nPc0I+i7gP4DdwIvA1UA3sB3YA2xLXpOKQS9J09OMoP9H4KfAauByYBDYRAj6i4FHk/NUatsVv/9+\n2k+SpLikXXWzGHgW+P0x1weBTwLDwFKgClw65jWTXnVTs3w5PPYYrFw5rbZKUuFlsepmBXAI+Bfg\n18A/AwuAHkLIkzz3pLwPYPlGkqajownvvxL4CvA08AAfLNOMJI8PGBgYOHtcqVSoVCoT3qwW9J//\n/LTbK0mFUq1WqVarqT4jbelmKfC/hJE9wLXA3YRSzqeAA0Av8DhNKN386EfhsWVLmiZLUnFlUbo5\nALxGmHQFuB7YBTwCbEyubQQeTnkfwLX0kjQdzdgCYR3wXWAOsA+4DWgHfgz8HvAK8GfAsTHvm/KI\n/vTpsAXC8DAsXJiy1ZJUQNMZ0Rdir5tGV10F3/42XHNNC1okSTlX2r1uGvmPkEjS1BQy6K3TS9Lk\nFS7onZCVpKkpXI3+8OHwzdhjx6Aty9ZLUgaiqNEvWRJW3Ozfn3VLJKkYChf04ISsJE1FIYPeOr0k\nTV4hg96VN5I0eQa9JJVc4VbdAJw6BYsWwZEj0NnZ5FZJUo5FseoGYPZsuOQS2LUr65ZIUv4VMujB\nCVlJmqzCBr11ekmaHINekkqukJOxEPakX706TMi6FYKkWEQzGQvQ0wNz5sDQUNYtkaR8K2zQg+Ub\nSZoMg16SSq7wQe/mZpI0scIHvSN6SZpYYVfdAJw8CV1dcPQozJvXpFZJUo5FteoGYO5cWLUKdu/O\nuiWSlF+FDnqwfCNJ51OKoHdCVpLGV/igd3MzSZpY4YO+NqJPOa8rSaVV+KDv7Q0hPzycdUskKZ8K\nH/RtbdbpJWkizQr6duBZ4JHkvBvYDuwBtgFdTbrPObnyRpLG16ygvwN4EahVyjcRgv5i4NHkvGWc\nkJWk8TUj6JcBfwp8l/q3tW4CNifHm4Gbm3CfcTmil6TxNSPovwl8FXi/4VoPUJseHU7OW2bNGtiz\nB957r5V3kaRi6kj5/huBg4T6fGWc14xQL+mMMjAwcPa4UqlQqYz3ERObPx/6+2FwMIzuJaksqtUq\n1Wo11Wek3dTs68CXgNPAPGARsAX4CCH4DwC9wOPApWPem3pTs0YbNsCNN8KttzbtIyUpd7LY1Owe\noA9YAWwAHiME/1ZgY/KajcDDKe9zXtbpJencmr2OvjZEvx+4gbC88rrkvKUMekk6t0LvR99o/364\n+mp4882mfaQk5U50+9E36uuDd9+Fgwezbokk5Utpgr62FcLOnVm3RJLypTRBD9bpJelcShf0bm4m\nSaOVLugd0UvSaKVZdQPw9ttwwQVw4gR0pP3OryTlUNSrbgAWLIBly8K+N5KkoFRBD5ZvJGmsUga9\nE7KSVFfKoHdEL0l1pQt6/7UpSRqtdEG/fDkcPw5vvZV1SyQpH0oX9LNmwdq1juolqaZ0QQ/W6SWp\nUSmD3jq9JNWVMugd0UtSXam2QKg5cQJ6e8Nze3tLbiFJmYh+C4SaRYugpwf27s26JZKUvVIGPVi+\nkaSa0ga9E7KSFJQ26B3RS1JQ6qB3czNJKumqG4AzZ2DxYhgaCs+SVAauumnQ3g4f/jDs3Jl1SyQp\nW6UNerBOL0lg0EtS6ZU+6J2QlRS70k7GAhw7Bn19YX/6WaX+I01SLLKYjO0DHgd2AS8AtyfXu4Ht\nwB5gG9CV8j7T0tUF3d3w8stZ3F2S8iFt0J8C7gQuA64B/gpYDWwiBP3FwKPJeSas00uKXdqgPwDs\nSI7/D9gNXATcBGxOrm8Gbk55n2mzTi8pds2sXPcDVwBPAj3AcHJ9ODnPhCN6SbHraNLnLAR+AtwB\n/HbMz0aSxwcMDAycPa5UKlQqlSY1p27dOrj33qZ/rCTNiGq1SrVaTfUZzVh1Mxv4L+BnwAPJtUGg\nQijt9BImbC8d876Wr7oBOH06bIEwPAwLF7b8dpLUUlmsumkDHgRepB7yAFuBjcnxRuDhlPeZto4O\nWL3arRAkxStt0H8cuBX4FPBs8lgP3A/cQFheeV1ynhnr9JJilrZG/z+M/4fF9Sk/u2n8R0gkxSyK\n74s6opcUs1JvgVBz+DCsXBm2RGjL8r9YklJyP/pxLFkSVty8+mrWLZGkmRdF0IPlG0nxiibonZCV\nFKtogt4RvaRYRRX0bm4mKUZRrLoBOHUKFi2CI0egs3PGbitJTeWqmwnMng2XXgq7dmXdEkmaWdEE\nPVinlxQng16SSi66oHdCVlJsopmMhbAn/erVYULWrRAkFZGTsefR0wNz5sDQUNYtkaSZE1XQg3V6\nSfGJMuit00uKSZRB74heUkyiC3o3N5MUm6hW3QCcPAldXXD0KMybN+O3l6RUXHUzCXPnwqpVsHt3\n1i2RpJkRXdCDE7KS4hJt0FunlxSLKIPeCVlJMYky6GulmwzmgiVpxkUZ9L29IeQPHMi6JZLUelEG\nfVubdXpJ8Ygy6ME6vaR4RBv0juglxaKVQb8eGAReAu5q4X2mxbX0kmLRqi0Q2oHfANcDQ8DTwBeB\nxu+jZrIFQs0770B3Nxw/Hvaol6QiyNMWCB8F9gKvAKeAfwc+16J7Tcv8+dDfD4ODWbdEklqrVUF/\nEfBaw/nrybVccUJWUgw6WvS5hfgq0pVXwu23w333Zd0SSWqdVgX9ENDXcN5HGNWPMjAwcPa4UqlQ\nqVRa1Jxzu/NO+MIXZvSWkjQlTzxR5cknq2fPv/WtqX9GqyZjOwiTsZ8G3gCeImeTsZJURNOZjG3V\niP408BXgF4QVOA8yOuQlSTMkun9hSpKKLE/LKyVJOWHQS1LJGfSSVHIGvSSVnEEvSSVn0EtSyRn0\nklRyBr0klZxBL0klZ9BLUskZ9JJUcga9JJWcQS9JJWfQS1LJGfSSVHIGvSSVnEEvSSVn0EtSyRn0\nklRyBr0klZxBL0klZ9BLUskZ9JJUcga9JJWcQS9JJWfQS1LJGfSSVHIGvSSVXJqg/3tgN/AcsAVY\n3PCzu4GXgEHgMynuIUlKKU3QbwMuA9YBewjhDrAGuCV5Xg98J+V9Sq9arWbdhNywL+rsizr7Ip00\nAbwdeD85fhJYlhx/DngIOAW8AuwFPpriPqXnL3GdfVFnX9TZF+k0a6T9F8BPk+MLgdcbfvY6cFGT\n7iNJmqKO8/x8O7D0HNfvAR5Jju8F3gN+OMHnjEy9aZKkZmhL+f4/B74MfBp4N7m2KXm+P3n+OXAf\nobzTaC+wMuX9JSk2+4BVM3Wz9cAuYMmY62uAHcAcYEXSqLR/oEiSMvAS8CrwbPL4TsPP7iGM2AeB\nP575pkmSJElqqfWE0f5LwF0Zt2WmfQ8YBnY2XOsmTHzvIXw/oSuDdmWhD3icUAJ8Abg9uR5jf8wj\nzGPtAF4E/ja5HmNf1LQTqgW1hR+x9sUrwPOEvngquZb7vmgnlHX6gdmEX+zVWTZohn0CuILRQf8N\n4GvJ8V3UJ7LLbinwB8nxQuA3hN+FWPujM3nuAJ4AriXevgD4G+AHwNbkPNa+eJkQ7I1y3xcfI6zE\nqdlEfaVOLPoZHfSDQE9yvDQ5j9HDwPXYH53A04RvnsfaF8uAXwKfoj6ij7UvXgZ+d8y1KfVFFlsT\nXAS81nDuF6rC/7Dh5HiY+v/AmPQT/qbzJPH2xyzC33CHqZe0Yu2LbwJfpf7te4i3L0YIf+g9Q1jO\nDlPsi/N9YaoV/PLUxEaIr48WAj8B7gB+O+ZnMfXH+4RS1mLgF4TRbKNY+uJG4CChJl0Z5zWx9AXA\nx4E3gQsIdfmxo/fz9kUWI/ohwiRcTR+jt0yI0TD1byD3En7JYzGbEPL/RijdQNz9AXAc+G/gD4mz\nL/4IuIlQsngIuI7w+xFjX0AIeYBDwH8S9g6bUl9kEfTPAB8i/FV9DmGny60TvSECW4GNyfFG6oFX\ndm3Ag4RVJg80XI+xP5ZQXzkxH7iBMKKNsS/uIQwAVwAbgMeALxFnX3QCv5McLyBs+76TgvTFnxBW\nWOylvr1xLB4C3iDsD/QacBthRv2X5HipVItcSyhX7KD+xbv1xNkfa4FfE/rieUJ9GuLsi0afpD4Q\njLEvVhB+J3YQliDX8jLGvpAkSZIkSZIkSZIkSZIkSZIkSZIk+H8YMQ5BvbqcBwAAAABJRU5ErkJg\ngg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0xb03a8bec>"
       ]
      }
     ],
     "prompt_number": 52
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print(\"Mass of propellant at t = 3.2s:\", results[2,t <= 3.2][-1])\n",
      "\n",
      "vmax = numpy.max(numpy.abs(results[1,:]))\n",
      "print(\"Maximum speed:\", vmax)\n",
      "print(\"Time at maximum speed:\", t[numpy.equal(numpy.abs(results[1,:]), vmax)])\n",
      "print(\"Height at maximum speed:\", results[0,numpy.equal(numpy.abs(results[1,:]), vmax)])\n",
      "\n",
      "hmax = numpy.max(results[0,:])\n",
      "print(\"Maximum height:\", hmax)\n",
      "print(\"Time at maximum height:\", t[numpy.equal(numpy.abs(results[0,:]), hmax)])\n",
      "\n",
      "timpact = t[numpy.less_equal(results[0,:], 0.)][2]\n",
      "print(\"Time of impact:\", timpact)\n",
      "print(\"Velocity at impact:\", results[1,numpy.equal(t, timpact)])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Mass of propellant at t = 3.2s: 36.0\n",
        "Maximum speed: 232.490174554\n",
        "Time at maximum speed: [ 5.001]\n",
        "Height at maximum speed: [ 537.75980961]\n",
        "Maximum height: 1350.44970068\n",
        "Time at maximum height: [ 15.787]\n",
        "Time of impact: 37.349\n",
        "Velocity at impact: [-85.98116603]\n"
       ]
      }
     ],
     "prompt_number": 69
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}