P3.4_Theis_pumping.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:51552bf5238a6b256f735732cfcc3a5980bde0a46621cd567cb3a9592b74ce62"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#Problem P3.4 Theis solution with two pumping\n",
      "\n",
      "In this notebook, we will work through one tutorial based the Theis solution for transient pumping, and investigate the superposition of drawdown from two interfering pumping wells. Two wells fully penetrate a 20-m-thick confined aquifer that is isotropic and homogeneous (Fig. P3.1). Storativity is estimated to be $2 x 10^{-5}$. The hydraulic conductivity is 100 m/d. The confining unit is composed of very low permeability material and is approximated as impermeable. Both wells have a radius of 0.5 m and are pumped continuously at a constant rate for 30 days; well A is pumped at 4000 L/min and well B is pumped at 12,000 L/min. Before pumping, the head is 100 m everywhere in the problem domain. The 800 m by 500 m problem domain in Fig. P3.1 is the near-field region of a problem domain that extends over many tens of square kilometers so that the aquifer effectively is of infinite extent and the composite cone of depression does not reach the boundaries after 30 days of pumping.\n",
      "\n",
      "We simpflied it to look like this:\n",
      "\n",
      "<img src=\"P3.4_figure.tiff\">\n",
      "\n",
      "Below is an iPython Notebook that builds a Theis function and plots results.  "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Problem 3.4, page 107 Anderson, Woessner and Hunt (2015)\n",
      "\n",
      "# import Python libraries/functionality for use in this notebook\n",
      "%matplotlib inline\n",
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "import scipy.special\n",
      "import sys, os\n",
      "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
      "\n",
      "# return current working directory\n",
      "os.getcwd()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 2,
       "text": [
        "'/Users/brclark/brc/tools/appliedGWmodeling/Chapter_3_problems-1'"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Define an function, class, and object for Theis Well analysis\n",
      "\n",
      "def well_function(u):\n",
      "    return scipy.special.exp1(u)\n",
      "\n",
      "def theis(Q, T, S, r, t):\n",
      "    u = r ** 2 * S / 4. / T / t\n",
      "    s = Q / 4. / np.pi / T * well_function(u)\n",
      "    return s\n",
      "\n",
      "class Well(object):\n",
      "    def __init__(self, x, y, rate, name):\n",
      "        self.x = float(x)\n",
      "        self.y = float(y)\n",
      "        self.rate = rate\n",
      "        self.name = name\n",
      "        self.swell = None\n",
      "        return"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Parameters needed to solve Theis\n",
      "r = 500      # m\n",
      "T = 2000     # m^2/d (100 m/d Kh x 20 m thick)\n",
      "S = 0.00002  # unitless\n",
      "t = 30.      # days\n",
      "\n",
      "#Q = pumping rate # m^3/d - but we'll enter it below in the well info"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Well information\n",
      "well_list =[]\n",
      "well_obj = Well(250, 250, 5760, \"Well A\")  # 4000 L/min = 5760 m^3/d\n",
      "well_list.append(well_obj)\n",
      "well_list.append(Well(550, 250, 17280, \"Well B\")) # 12000 L/min = 17280 m^3/d\n",
      "\n",
      "# Grid information as requested in problem\n",
      "x = np.linspace(0, 800., 50)  # x-direction 0 to 800 m, 50 m increments\n",
      "y = np.linspace(0, 500., 50)  # y-direction 0 to 500 m, 50 m increments\n",
      "xgrid, ygrid = np.meshgrid(x, y) # make a grid with these coordinates"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "###We want to explore drawdown as a function of time\n",
      "So, set up an array of times to evaluate, and loop over them. Also, we can specify a distance from each well at which to calculate the curve of drawdown over time."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "times = np.linspace(0.,30.,31) # linear interpolation of time from 0 to 30 days, make 30 increments days at 0.5\n",
      "rdist = 25  # this sets the distance to plot drawdown over time\n",
      "print times"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[  0.   1.   2.   3.   4.   5.   6.   7.   8.   9.  10.  11.  12.  13.  14.\n",
        "  15.  16.  17.  18.  19.  20.  21.  22.  23.  24.  25.  26.  27.  28.  29.\n",
        "  30.]\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "###We will want to normalize our plots\n",
      "Let's figure out the maximum drawdown to use for setting our colorbar on the plots."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#let's find the maximum drawdown\n",
      "drawdown_grid_max = np.zeros(xgrid.shape, dtype=float)\n",
      "for well_obj in well_list:\n",
      "    r = ((well_obj.x - xgrid)**2 + (well_obj.y - ygrid) ** 2) ** 0.5\n",
      "    s_max = theis(well_obj.rate, T, S, r, times[-1])\n",
      "    drawdown_grid_max += s_max\n",
      "max_drawdown = np.max(drawdown_grid_max)\n",
      "print max_drawdown"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "15.4068299862\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "###Clobber the PNG output files\n",
      "We will make a sequence of PNG files from which to render the movie. It's a good idea to delete the existing ones from our folder. \n",
      "\n",
      "**do this kind of thing with caution!** \n",
      "\n",
      "**You are deleting files and will _not_ be asked to confirm anything!**"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for cf in os.listdir(os.getcwd()):\n",
      "    if cf.endswith('.png'):\n",
      "        os.remove(cf)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "###Loop over time and make figures\n",
      "We will make a figure with the drawdown contours over the whole grid, and another with the addition of a plot of drawdown over time at a distance of `rdist` from each well. NB --> this second calculation doesn't take account of both wells. We could do so, how might we?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Note that this section of code is saving figures for animation - not plotting them!\n",
      "\n",
      "from IPython.display import clear_output\n",
      "# to make our plots of drawdown over time a one point, we can\n",
      "# predefine the response as np.nan. That way, when we plot incrementally\n",
      "# as we calculate through time, only the times for which calculations\n",
      "# have been made will appear using plt.plot()\n",
      "for well_obj in well_list:\n",
      "    well_obj.swell = np.ones_like(times)*np.nan\n",
      "    \n",
      "\n",
      "# using \"enumerate\" we get both the iterant (t) and a counter (i)    \n",
      "for i,t in enumerate(times):\n",
      "    # the following stuff just writes out a status message to the screen\n",
      "    clear_output()\n",
      "    perc_done = (i/float(len(times)-1)) * 100\n",
      "    sys.stdout.write('working on time {0}: {1:2.2f}% complete'.format(t,\n",
      "                                                                      perc_done))\n",
      "    if i < len(times):\n",
      "        sys.stdout.flush()\n",
      "    # here's the end of the silly shenanigans of plotting out status to the screen    \n",
      "\n",
      "    # now we calculate the drawdown for each time.\n",
      "    drawdown_grid = np.zeros(xgrid.shape, dtype=float)\n",
      "    for well_obj in well_list:\n",
      "        r = ((well_obj.x - xgrid)**2 + (well_obj.y - ygrid) ** 2) ** 0.5\n",
      "        s = theis(well_obj.rate, T, S, r, t)\n",
      "        well_obj.swell[i] = (theis(well_obj.rate, T, S, rdist, t))\n",
      "        drawdown_grid += s\n",
      "     \n",
      "   \n",
      "    # drawdown field and at a point \n",
      "    plt.subplot(1, 3, 1, aspect='equal')\n",
      "    im = plt.contourf(xgrid, \n",
      "                 ygrid, \n",
      "                 drawdown_grid, \n",
      "                 np.linspace(0,max_drawdown,10))\n",
      "    # crazy stuff to make the colorbar fit better. very optional!\n",
      "    divider = make_axes_locatable(plt.gca())\n",
      "    cax = divider.append_axes(\"right\", \"5%\", pad=\"3%\")\n",
      "    plt.colorbar(im, cax=cax).ax.invert_yaxis()\n",
      "    for well_obj in well_list:\n",
      "        plt.text(well_obj.x, well_obj.y, well_obj.name)\n",
      "    plt.title('Drawdown at time = {0:.0f}'.format(t))\n",
      "    \n",
      "    # also plot drawdown at a point from Well A\n",
      "    ax = plt.subplot(1,3,2)\n",
      "    plt.hold(True)\n",
      "    plt.plot(times,well_list[0].swell)      #remember 0 index in Python\n",
      "    plt.plot(t,well_list[0].swell[i],'ro')\n",
      "    ax.set_xlim([0,times[-1]])\n",
      "    ax.set_ylim([0,max_drawdown])\n",
      "    plt.tight_layout()\n",
      "    plt.savefig('s{0}.png'.format(i))\n",
      "    plt.close()\n",
      "    \n",
      "    # also plot drawdown at a point\n",
      "    ax = plt.subplot(1,3,3)\n",
      "    plt.hold(True)\n",
      "    plt.plot(times,well_list[1].swell)\n",
      "    plt.plot(t,well_list[1].swell[i],'bo')\n",
      "    ax.set_xlim([0,times[-1]])\n",
      "    ax.set_ylim([0,max_drawdown])\n",
      "    plt.tight_layout()\n",
      "    plt.savefig('spt{0}.png'.format(i))\n",
      "    plt.close()\n",
      "    \n",
      "# Let's finish with a drawdown only plot --> make a second set of figures with only the \n",
      "    # make a plot    \n",
      "    plt.subplot(1, 1, 1, aspect='equal')\n",
      "    im = plt.contourf(xgrid, \n",
      "                 ygrid, \n",
      "                 drawdown_grid, \n",
      "                 np.linspace(0,max_drawdown,10))\n",
      "    plt.colorbar().ax.invert_yaxis()\n",
      "    for well_obj in well_list:\n",
      "        plt.text(well_obj.x, well_obj.y, well_obj.name)\n",
      "    plt.title('Drawdown at time = {0:.0f}'.format(t))\n",
      "    plt.savefig('s_only{0}.png'.format(i))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "working on time 30.0: 100.00% complete"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD3CAYAAAA0Vx7KAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcXFWd9/HPlwCGkJAIYYBAIJFhSUAhjERkGRsHIS6A\nOiLBQREchhnWB32UxDU8jriMjuiADqPARGQVhUlE2ZRAECSgCYYkDSIECEtYEzbjBPJ7/rinkkql\nqrqqu7rqVtX3/XrVq+9ybt1fd93761Pn3nuOIgIzM8uvjVodgJmZVedEbWaWc07UZmY550RtZpZz\nTtRmZjnnRG1mlnNO1GZmfZA0RVKvpD9KOqtCme+m9fdKmlTPtn1xojYzq0LSEOA8YAowEThG0oSS\nMu8B/joidgH+Cfh+rdvWwonazKy6ycCDEbE0IlYDVwBHlpQ5ApgJEBF3AaMkbVvjtn1yojYzq257\n4LGi+WVpWS1lxtSwbZ+cqM3Mqqu1nw0NVgAbD9Ybm5nlhaS6OjWKiOKk+zgwtmh+LFnNmCpldkhl\nNqlh2z45UZtZl7ilxnIHly64B9hF0jjgCeBo4JiSMrOAU4ErJO0HrIiI5ZKeq2HbPjlRm5lVERGv\nSToVuAEYAlwYEUsknZTWXxARv5D0HkkPAq8Ax1fbtt4Y5G5OzazTZU0ftdeoS5o+Ws4XE83Mcs6J\n2sws55yozcxyzonazCznnKjNzHLOidrMLOecqM3Mcs6J2sws55yozcxyzom6y0laI+lNrY6jUSTt\nKOklSbl6ssxsIJyom0zSUkmvSnpR0guSfiPpJCeWDUmaIemSPsoslfTOwnxEPBoRIyLHfSNI+rGk\nJ9Mx8JCkz5Ws/7s0dNMrkn4tacdWxWr54ETdfAG8LyK2AHYEvgacBVxYaQNJ/pwqCwaxH+BB8lVg\nfDoG3g2cJmkKgKTRwE+BzwFvJOu57cpWBWr54ATQQhHxUkTMJuv68DhJEwEk/bek70v6haSXgR5J\n75U0X9JKSY9K+lLhfSTNlPTJNL19as44Oc3vnLpaLJT9tKQnJC2TdEJxPJJGSvqRpKdTTfVzhZq+\npEck7ZOm/yHtY0Ka/4Ska9L0DElXpZhelHSfpL+p9DeQ9J30+6yUdI+kA9PyKcB04OjUlDG/zLaX\nkP2zm53K/F9J41JsG6UycyR9OX1zeUnSLEmjJV2a9jlP0k5F77m7pJskPZdqtUfV/IHWKCIWRcSq\nokWvAU+n6Q8C90XETyPif4EZwF6Sdm10HNY+nKhzICLuJutM/KCixccAX46I4cBvgJeBYyNiJPBe\n4F8kFcZemwP0pOl3AA8Bf1s0fxusTX6fAg4Bdk0/i/0HMAIYn7b7GKm7xjL7+FP6WZifU/Q+hwOX\nAyPJ+uk9r8qvPw/Yi6z2eBnwE0mbRsT1wDnAFakpY1LphhHxUeBRsm8oIyLimxX2cTRwLNkQSDsD\nd5J9g9kSWAJ8CUDS5sBNwI+BrYGpwPcqDUYq6Xup+arca0GV37mw7SvAIuBfI+L3adUewL1Fv+Or\nwIPAntXezzqbE3V+PEGWOAqujYg7ASLiLxFxa0QsSvMLyQbJLCTK24ADU+33IOAbwAFp3TuAW9P0\nh4GLImJxSgDFtfIhZAltekS8EhGPAN8CPpqK3Fq0vwPJvr4X5v+2aB8AcyPi+tRO/GOyRFxWRFwa\nES9ExJqI+HfgDcBuhbAYeLNGABdHxMMR8SLwS+CBiPh1RLwO/AQo/BN4H/BwRMxM8SwAfgaUrVVH\nxMkR8cYKr72rBhVxMjCc7J/lv0qanFZtDrxYUvzFVNa6lBN1fuwAPJ+mg/UHxETS2yTdkpolVgAn\nAVsBRMSfyDor35ssUf8ceCJ9XS5OotuVvO+jRdOjyYYNeqRkfWEgztuAg5SNrDyELMEdkJoNRqak\nVrC8aPpVYGildvbUXLFY0gpJL5DVwkeXKzsAxfGsYl0zQ2G+kAR3At5WXDMGPgJs0+B4AIjMHLK/\nZWHUj5eBLUqKjgReGowYrD04UeeApH3JRiu+vUqxy4BrgR0iYhTwn6z/+d1KVvPbJCKeSPMfJ2tS\nKCTRJ8nadAuKp58FVgPjStYvA4iIB8mS7mnArRHxEvAU8E/A3KJtar7bQtJBwKeBoyJiVES8EVjJ\nulp0Le9V790d1co/Sva7FdeMR0TEKeUKS/rP1O5d7rWwjpg2IftHC1lTyNpvIKk5Zue03FpA0kWS\nllf6TFNlY356LZT0mqRRtWxbKyfq1ihcoNtC0vvI2nMvKTRtUP7r/nDghYj43/Q1+SOsn3RuJRuz\n7bY0PyfNzy26Ve0q4OOSJkgaRlHTR2oGuAr4iqThqaZ8JlnTRek+CjX0OSXzlWKvZATZhbRnJW0q\n6YusX5t8ChhXuKBZwXKyRFaNKkyXug7YVdKxkjZJr30l7V6ucET8c0rk5V5vLhuItLWkqZI2lzRE\n0mFk/2D/JxW5BthT0gclDSX7jBZExAN9/I42eC4GplRaGRHfjIhJ6TrKdGBORKyoZdtaOVG3xmxJ\nL5LV4KaTtQUfX7Q+2LDmdzLw/9J2X2DDW7ZuI0vmhUT9G2CzonnSBbpzgV8DDwC/KtnPaWQ1u4fI\nasmXkh1oBbeW7KN0vlLslWqx16fXA8BS4M+s3xzzk/TzOUn3VHiPrwKfT00Vn6ywvyiZLrs+fUs4\nlOwi4uNk30C+CmxaYd/9EcA/k31TeQ74MvDRdEGZiHgW+HvgK2RNYW9N8ViLRMRc4IUai3+ErOLV\nn20rqmnMRElLyS5ovA6sjojJkrYkSxY7kZ1kHy78F5E0HTghlT89Im4caKBmZv2lAY6ZqGwU8dmV\nvimlMsPIrgHtXFSjrmnbvtQ6CnkAPRHxfNGyacBNEfENSWel+WnK7gU+GphIdiHqZkm7RsSa/gZp\nZjZgE3rKL39lDrw6Z938s/3ew+HA7cVJulFqTdSwYdveEay7PWsmWXvlNOBI4PKIWA0sVTZ8+mTg\ntwML1cxsEGzek70Knj27v+80laJmj0aqtY06yGrG90g6MS3bJiIKtz0tZ90tTGNIdwoky1h3i5eZ\nWceRNJLsVtj/6atsf9Raoz4gIp6UtDVwk6Te4pUREVkbUEXrreujrJnZekrbjJtJ0uVkrQejJT1G\ndifOJimuC1Kx9wM3RMSfK2y7Vdr2ixFRfIG+JjUl6oh4Mv18RlmfDpOB5ZK2jYinJG3HuocIHgfG\nFm2+Q1q2niWxU+mipjtvxgpOnTGq62NoRBy7z36k70J9mHEZzPjIgN+mI+IYSAy9hzfm3MrLsTlB\nAz+2BiIijqmhzEyyJuC6t61Fn00fkoZJGpGmNye7fWkhWR8Ox6Vix5E9jEFaPjXdFzse2IWsPwfr\nQLvPfqQhSdoax59J56mlRr0NcE165mBj4NKIuDHd13qVpE+Qbs8DiIjFkq4CFpM9zHBynvsGtvo5\nCbSHwufUqBq2tU6fiToiHibrQ6J0+fNs2PtaYd05ZD2f5drknqGtDiEXMUBtcQx2gu7p912mjZWH\nOBoZw+6zH+lXss7LsWk1PvDS8J1KkYc2aquNa9Cdox1r1xP0yIAvJkoKJtSY65aopRcvy6nnPmrr\nMk7QncfNIe3JfX3YBnwxqvP5820vrlHbWj55u4tr1+3DNWoDnKS7mT/7/HONusv5JDVw7TrvXKPu\nUm6HtnJ8TOSTE3WXcYK2vvj4yB8n6i7iE9Bq5WMlX9xG3QV80ll/uN06P1yj7nBO0jZQPoZaz4m6\nQ7kt2hrJx1JrOVF3IJ9UNhh8XLWOE3UHcS3aBls3Hl+Shkq6S9ICSYslfbVK2X0lvSbpg2l+rKRb\nJC2SdJ+k0/sTgxN1h+jGE8hao9uOtYhYBRwcEXsDbwEOlnRgaTlJQ4CvA9ezbjDw1cCZEbEHsB9w\niqQJ9cbgRN3mXIu2Vui2Yy4iXk2TmwJDgOfLFDsNuBp4pmi7pyJiQZp+GVhCNgB4XZyo21i3nSyW\nL910/EnaSNICYDlwS0QsLlm/PXAk8P20aIPOryWNAyYBd9W7f99H3aa66SSx/Orv6DEt8X8rLO+d\nA/fPWTe/ZMMiEbEG2FvSSOAGST0RUbQR5wLTIiKUjVu43sADkoaT1bbPSDXruniElzbjBG15NJjJ\numEjvFxYY677RPURXiR9AfhzRHyzaNlDrEvOo4FXgRMjYpakTYCfA7+MiHP7E7+bPtqIk7TlVScf\nm5JGSxqVpjcD3gXMLy4TEW+KiPERMZ6s5vwvKUkLuBBY3N8kDU7UbaOTTwTrDB18jG4H/Dq1Ud8F\nzI6IX0k6SdJJfWx7AHAs2Z0i89NrSr0BuI26DXTwCWCWexGxENinzPILKpQ/vmj6dhpQIXaNOuec\npK2d+HgdHE7UOeX7o61d+bhtPCfqHPKBbu3Ox3Bj1ZSoJQ1JjeCz0/yWkm6S9ICkGwtXRNO66ZL+\nKKlX0qGDFXin8gFuncLHcuPUWqM+A1jMuqdtpgE3RcSuwK/SPJImAkcDE4EpwPckudZeIx/YZlZO\nn0lU0g7Ae4Afsu6G7iOAmWl6JvD+NH0kcHlErI6IpcCDwORGBtypnKStE/m4boxaarvfBj4NrCla\ntk1ELE/Ty4Ft0vQYYFlRuWXA9gMNstP5YLZO5uN74KreRy3pfcDTETFfUk+5MunZ9mrPZpZdd96M\nFWunJ/cMZXLP0L6j7UA+iM3WN2/OKubNWdXqMHKlrwde9geOkPQeYCiwhaRLgOWSto2IpyRtBzyd\nyj8OjC3afoe0bAOnzhhVbnFXcZK2blFP502lFbfzz145WGG1japNHxHx2YgYm55fnwr8OiI+CswC\njkvFjgOuTdOzgKmSNpU0HtgFmDc4obc3J2nrNj7m+6/eR8gLzRhfA66S9AlgKfBhgIhYLOkqsjtE\nXgNOjlZ0z5dzPmDNrB7u5rTJnKSt29XbJWreujltBd/j3ERO0mbWH07UZtZUrrDUz4m6SXxwmll/\nOVE3gZO02fp8TtTHiXqQ+YA0a3+Slkr6Q+qcboNbjiXtLulOSaskfarM+vU6tquXE/UgcpI2q6zN\nzo8AeiJiUkSU67/oOeA04Jtl1sGGHdvVxYnazKw2FW/Zi4hnIuIeYPUGG5Xv2K4uTtSDpM1qC2Yt\n0UbnSQA3S7pH0ol1bluuY7u6eHDbQdBGB59Z19jphN6yy1fNmceqOeuanSv0LHJARDwpaWvgJkm9\nETG3r33W0rFdLZyozayl6umwaTAM7ZnM0J51zc4rzz5/gzIR8WT6+Yyka8j62e8zUVO+Y7sfRcTH\n6onRTR8N5tq0WWeRNEzSiDS9OXAosLBS8eKZCh3b1ZWkwYm6oZykzfon5+fONsBcSQuAu4CfR8SN\nkk6SdBKApG0lPQacCXxe0qOShpd5r37d9eGmDzOzKiLiYWDvMssvKJp+ivX74i/3PrcCt/YnBteo\nGyTnNQIza2NO1GaWC67sVOZE3QA+wMxsMLmNeoCcpLvIhndt1e6UhkVhXciJ2qyagSTnSu/jpF1R\nq++pzisn6gFwbbpDNSo51/L+TtpWAydqs4LBTtDV9umEbVX4YqIZtCZJ52n/OeJvqhtyjbqffDB1\niDwlSNeurQLXqK175SlJF8trXNYyTtT94Np0B8h7Msx7fNZUVRO1pKGS7pK0QNJiSV9Ny7eUdJOk\nByTdKGlU0TbTJf1RUq+kQwf7FzCry/m0TxJslzgHgStD66uaqCNiFXBwROwNvAU4WNKBwDTgpojY\nFfhVmkfSROBoYCIwBfieJNfaLR/aMfG1Y8zWcH0m0Yh4NU1uCgwBXgCOAGam5TOB96fpI4HLI2J1\nRCwFHiTrYLtj+D99m2rnhNfOsVtD9JmoJW2U+mFdDtwSEYuAbSJieSqynKy/VoAxwLKizZcB2zcw\nXrPu5GTd1fq8PS8i1gB7SxoJ3CDp4JL1IalaZ9hl1503Y8Xa6ck9Q5ncM7S2iM3q1SlJ7ny64ta9\nOQuz17O/W9F34SZJ1+F+COxBltNOiIjfFq0fDfwY2JYsr34zIv47rVsKvAi8DqyOiLpbGWq+jzoi\nVkq6DvgbYLmkbSPiKUnbAU+nYo+zfufZO6RlGzh1xqhyi3PNzR5tqFOSdBfpeXP2gpX0Hr4T559d\nYbjZ5voO8IuI+JCkjYHNS9afCsyPiOkpad8v6ccR8RpZYu+JiOf7u/O+7voYXbijQ9JmwLuA+cAs\n4LhU7Djg2jQ9C5gqaVNJ44FdgHmYtUInJulO/J1yLrUmHBQRFwFExGsRUfrf40lgizS9BfBcStJr\n32YgMfRVo94OmJnu3NgIuCQifiVpPnCVpE8AS4EPp19gsaSrgMXAa8DJEdGvMcLMrIIuaQLJkfHA\nM5IuBvYCfgecUXSjBcAPgF9LegIYQcqJSQA3S3oduCAiflBvAFUTdUQsBPYps/x54JAK25wDnFNv\nIGYN5ZqnlfgYl9RU7ssbLtqYLA+eGhF3SzqX7JbkLxaV+SywICJ6JO0M3CRpr4h4CTggIp6UtHVa\n3hsRc+uJ3X191Mjt05YrrlU3zNI5j/DInKrn9zJgWUTcneavJj07UmR/4CsAEfEnSQ8DuwH3RMST\nafkzkq4hu2XZidq6nGvTVodxPTsxrmfdYAW3nX37euvTTROPSdo1Ih4ga01YVPI2vWn5byRtQ5ak\nH5I0DBgSES9J2hw4FDi73hidqM3M+nYacKmkTYE/ASdIOgkgIi4ga+69WNK9ZNfzPhMRz0t6E/Az\nSZDl20sj4sZ6d+5EbTYI5twAPYcN8k66oPkjL02OEXEvsG/J4guK1j8LHF5mu4eAvQe6fydq6yxN\nbvaYc0Pf6wY9YVvHc4dJHeDMH8J3Zq2bP+xLcOJ56+Y/dSF8+38qb//xc+Gnd2TTPZ+F3z1Yvtyz\nL8ImH4ALrh94zJ2gUpI+n+xqU8G+N8CJ962b/1QvfHtp5ff9+EL46VPZdM88+F2Z5z165sHut8Ok\nM2DiKfCDKv8wrP05UXeAAyfCHb3Z9Jo18NxLsPjRdevvvB8OmFB5e2nd3fhS9irnJ7fDlH3g8tsa\nEnZbq1aT3pN1V5rWkD07fOfj67a5cyUcUOXBXLHuMyieLi1z2Vtg/nfgN1+Hs2bCa6/X+UtY23Ci\nrkFe2skqeftucGdK1IsehT13hBGbwYqX4S+rYcljsM/OWU2557Pw1k/ClC/BUy/Ut58r5sK/HgtP\nr4THn2v879EuqiVpyDqDWJyml5I9LTEMeBm48QZY8jLss0VWU+6ZB2+9E6bcA0/9pb44IoDz4cVX\nYfhQGOKzuWO5jboDjNkKNh4Cjz2TJey3754l0jvvhy02g7eMy8qd9l8w+/Ow1RZw5Vz43CVw4em1\n7eOxZ7IEvdd4+ND+2faffH/f2zVVTm7LG03WH/DTZDXricCzaXoY8JYRWbnTemH2JNhqU7jySfjc\nH+HCPWvbRwD/sBDeIPjjLfCdEyt/E7L250TdIfbfPWv+uKMXPnkkPP483LEERm6eNXvc/3hW2z7k\nC1n519fAmC1rf/8rb88SNMBRB8AJ381hom6CvmrTBXsA95El56NYl6g3B3Z4Ae5/FRa9DIfck5V/\nPWDMG2qPo9D0sc8W8OxHYf/PwGH7wI5b1/4e1j6cqDvEARPgN0tg4SPw5nEwdmv45jUwchic8K7s\na/IeO8Id3+jf+19+GyxfAT+ek80/+QI8+AT89ZhG/QadZU+yRP0Q8Cbgr4ArgeHAu0mfx3C4420D\n39foLbKmrbvud6LuVG7V6hD7T4Cf3w1bjci+Ar9xOKx4JWv+2H932HUMPLMSfpvasle/tv4Fx2oe\neBxeWQXLLoaHf5i9pv29LypWswfwW2AkWe13BFkb9aK0btfN4Zn/hd+mLpdXr4HFL9e3j0J3Z6/+\nBeY/BH+9XWNit/xxjbpD7LljdrfHsT3rlr1lXHYSb5naRK+eBqf/F6x8NbtD4MwjYOKOfb/3FXPh\ng29ff9nf7w9T/w2+MLVRv0FnGU92t0dxz2U7A38h6wNz043g6r3g9F5Y+Rq8FnDmTjBxeO37+IeF\nsNlG8JfFcPzfwaSdG/kbWJ6oFb2QSoolsVPfBXMi73d9WNKEi4m1tlH3paEPwXT404k6AiJiQJdK\nJcUX4rM1lf2yzhnw/hrNTR9mZjnnRG1Wh0bUhP1IudXLidrMLOecqGvQe3j7tKd3tSa11Q6kRuza\ntPWHE7VZP/Qn4Q5Kku7wC4mWcaI266d6Eq9r0jYQTtRmA9BzWPUk3Nd6aw+SpktaJGmhpMskvaFk\nfY+klZLmp9fni9aNknS1pCWSFkvar979+4EX6yyn0JLOmVqSjN3s0RSSxgEnAhMi4i+SrgSmAjNL\nit4aEUeUeYvvAL+IiA9J2pisy5e6uEZtZrmVkwv5LwKrgWEp0Q4DHi9TboOHZCSNBA6KiIsAIuK1\niCgzFER1TtTWeVzTtAaKiOeBbwGPAk8AKyLi5tJiwP6S7pX0C0kT0/LxwDOSLpb0e0k/SCOT18WP\nkNfBj5K3kZz0TT1ouuSfUe/hOzFBjzTkEfJac07p/iTtDMwGDgJWAj8Bro6IS4vKjABej4hXJb0b\n+E5E7CrprcCdwP4Rcbekc4EXI+KL9cTfZxu1pLHAj8h6agzgvyLiu5K2JOu5cSeygSw+HBEr0jbT\ngROA14HT+zM8utmAtKit2trPvDmrmDdnVbUibwXuiIjnACT9DNgfWJuoI+KloulfSvpeypHLgGUR\ncXdafTUwrd4Ya7mYuBo4MyIWSBoO/E7STcDxwE0R8Q1JZ6WdT0tV/qPJBrbYHrhZ0q4Rsabe4Mys\njC6pTTfL5J6hTO4Zunb+/LM3aELuBb4gaTNgFVmniPOKC0jaBng6IkLSZLLWiufTusdSDnwgbbuI\nOvXZRh0RT0XEgjT9MrCELAEfwbqrnjOBwngfRwKXR8TqiFgKPAhMrjcwswFzQrMGiIh7yVoV7gH+\nkBb/QNJJkk5K8x8CFkpaAJxLdldIwWnApZLuBd4CnFNvDHXdnpduU5kE3AVsExHL06rlwDZpegxZ\nn+kFy8gSe9vrPXwnt1Nba/mfT0tExDeA0vGRLihafz4VGttSot93IPuvOVGnZo+fAmdExEsqGkkz\nVferXZXcYN15M1asnS796mHWMJ3UVt0lSXrOwuz17G4j4Xcr+t6gC9SUqCVtQpakL4mIa9Pi5ZK2\njYinJG1HNugyZPcXji3afAfK3HN46oxR/Y/arB6dkKy7JEkD9Lw5e/UenuWIMm3GXafPNmplVecL\ngcURcW7RqlnAcWn6OODaouVTJW0qaTywCyUN72ZN10WJzjpPLQ+8HAAcCxxc9Bz7FOBrwLskPQC8\nM80TEYuBq4DFwC+Bk6MVN2sPkpw8KWX90a7Jul3jtobps+kjIm6nckI/pNzCiDiHflzZNBt07dYM\n4iRt+BHyfnGtus21S/JrlzgHgc+x9TlRW3fKcxI8hXzHZ03nRG3dK48JMW/xWC44UfeTv5p1kLwk\n7DzEkAM+tzbkgQPMClpxodHJ2WrgRD0AfqS8AxUnzsFM2k7QVgcnarNKGp20nZytn5yoB8i16i5R\nKclWSuBOyv3i9unynKjNBsIJ2ZrAd300gGsBZgPn86gyJ2ozs5xzom4Q1wbMOpOkiyQtl7SwaNmW\nkm6S9ICkGyVt0G+zpLGSbpG0SNJ9kk4vWndFUSd3D0uaXy0GJ+oGcrI265+cnzsXA1NKlk0jGzN2\nV+BXlB+wtjDe7B7AfsApkiYARMTUiJgUEZPI+vr/abUAnKjNzKqIiLnACyWLK40ZW7xdufFmxxSX\nSf39fxi4vFoMTtQNlvOagZk1RqUxY8sqGW+22EHA8oj4U7XtfXveIPC91Wa1a1blptI5WRijsb/6\nGjM2jTd7Ndl4sy+XrD4GuKyvfThRm1lXK4zRWHD2FTVtVmnM2PUUjTf746LxZgvrNgY+AOzT187c\n9DFI3ARi1rc2Pk8qjRm7VpXxZgsOAZZExBN97cyJehC18UFoZomky4E7gN0kPSbpeCqMGStpjKTr\n0qblxpt9d9FbH00fFxHXxtCKcWclxZLoniTm9mqzDdVakZmgR4gIDWRfkiJm1Vj2CAa8v0ZzjdrM\nLOecqJvATSBm6/M5UR8n6ibxgWlm/eVE3URO1mY+D/qjz0Rdb4ckkqZL+qOkXkmHDlbg7coHqXUz\nH//9U0uNuuYOSSRNJLvlZGLa5nuSXGsv4YPVzOrRZxKts0OSI4HLI2J1RCwFHgQmNybUzuJkbd3G\nx3z/9be2W6lDkjHAsqJyy4Dt+7mPjucD18xqMeC+PvrqkAQou+68GSvWTk/uGcrknqEDDaUtuQMn\n6wb1VErmzVnFvDmrBjGa9tPfRF2pQ5LHgbFF5XZIyzZw6owNBkToWk7W1snq/eZYWnE7/+yVjQ6p\n7fS36aNShySzgKmSNpU0HtgFmDewELuDm0GsE/m4boxabs+ruUOSiFgMXAUsBn4JnByt6EykTfmg\nNrNy3ClTDrkZxDpBoyoe7pTJTybmkmvW1u58DDeWE3VO+UA3swIn6hxzsrZ25OO28Zyoc84HvbWT\nTjxeK/R3dJSkRZJel1RxzENJZ0haKOk+SWcULZ8saV4a9eVuSftWi8GJug30Hr5TR54A1lk6+Bgt\n19/RQrKBaW+rtJGkPYF/BPYF9gLeJ2nntPobwBciYhLwxTRfkRN1G+ngE8HaXCcfm+X6O4qI3oh4\noI9NdwfuiohVEfE6cCvwwbTuSWBkmh5FhQcDC5yo20wnnxDWnnxMVnQfcFDqFnoY8F6yp7Uh63H0\nW5IeBf4NmF7tjQbc14c1nx85t7xoqyR9fvnFc57PXo0WEb2Svg7cCLwCzAdeT6svBE6PiGskHQVc\nBLyr0nv5gZc254RtrdKsJN2wB14Oq7HsDRs+8CJpHDA7It5csvwW4FMR8fsaYjgHeDQi/lPSixGx\nRVouYEVEjKy0rZs+2lxb1WisY/i4W0/FfyKS/ir93JHs4uNladWDkt6Rpt8JVG3vdtNHByicNK5d\nWzN0W5JO/R29Axgt6THgS8DzwH8Ao4HrJM2PiHdLGgP8ICLemza/WtJWwGqyvo9eTMv/CThf0huA\nP6f5yjFiXAg8AAAIU0lEQVS46aOzOFnbYGpFks5D00eruemjw/ieaxssPq5ax4m6Q/mkskby8dRa\nTtQdzLVrawQfQ63nRN0FfKJZf/nYyQcn6i7h2rXVy8dLfjhRdxmffNYX/1PPHyfqLuQT0SrxcZFP\nTtRdzAnbivlYyC8/mWh+srHLOUHnn2vUtpZP2O7jz7w9uEZt63Htujs4QbcX16itLLdfdy5/ru1n\nUGrUkqYA5wJDgB9GxNcHYz82+FzD7hxO0O2r4TVqSUOA88gGg5wIHCNpQqP30wjz5qxqdQi5iAH6\njqNQwx7Mk33Owr7LNEMe4mhkDP393PJybNrgNH1MBh6MiKURsRq4AjhyEPYzYHk4EPMQA9QXx2Al\n7DwkSMhHHI2KYSCfU16OTRucpo/tgceK5pcBbxuE/ViLuVkkv9zM0VkGo0bd/JEIrKV84TE//Fk0\nnqSLJC2XtLBo2Zcl3StpgaRfSRpbYdtRkq6WtETSYkn7peV7SbpT0h8kzZI0omoMjR7hJQUyIyKm\npPnpwJriC4qSnMzNrGatHOFF0kHAy8CPCoPbShoRES+l6dOAvSLiH8vsdyZwa0RcJGljYPOIWCnp\nbuCTETFX0vHA+Ij4YqWYBqPp4x5glzRq7xPA0cAxxQXyNsyNmVklKZmOK1n2UtHscODZ0u0kjQQO\niojj0javASvT6l0iYm6avhm4HqiYqBve9JGCORW4AVgMXBkRSxq9HzOzVpL0FUmPAscBXytTZDzw\njKSLJf1e0g8kDUvrFkkq3GRxFFC26WTtvloxuK2ZWTNJilsqrFuQXgUz2fBbf6pRzy40fZSsmwbs\nFhHHlyx/K3AnsH9E3C3pXODFiPiipN2A7wJbAbOA0yNidKX4/Qi5mXW1vdOrYGb9b3EZ8Isyy5cB\nyyLi7jR/NTANICLuBw4DkLQr8N5qO2j6I+SSpkjqlfRHSWcN4n7KXandUtJNkh6QdKOkUUXrpqeY\neiUd2qAYxkq6RdIiSfdJOr1FcQyVdFe6Qr1Y0ldbEUd63yGS5kua3cIYlqar7fMlzWtFHGXuBnhb\nC2LYLf0NCq+Vkk5vQRzT0zmyUNJlkt7QiuOizph3KZo9EphfWiYingIeS4kY4BBgUdp+6/RzI+Dz\nwPer7jAimvYie6T8QWAcsAnZN44Jg7Svg4BJwMKiZd8APpOmzwK+lqYnplg2SbE9CGzUgBi2BfZO\n08OB+4EJzY4jvfew9HNj4LfAgS2K45PApcCsVnwm6b0fBrYsWdbsY2MmcELRZzKyFX+Long2Ap4k\nayttWhzpfR4C3pDmryRr821oDEDcUuMLiJJtLye7MeJ/yZ4ROYGsdrwwxfJT4K9S2THAdUXb7gXc\nDdwL/AwYmZafTpYP7gfO6TP+Rn7YNfyx3g5cXzQ/DZg2iPsbx/qJuhfYJk1vC/Sm6enAWUXlrgf2\nG4R4riX7r9qyOIBh6cDZo9lxADuQXeE+mKy9ryWfCVmi3qpkWdPiIEvKD5VZ3srj4lBgbgv+Flum\nZPVGsn9Ys4F3NToGBpCo8/BqdtNHuacWt2/i/reJiOVpejmwTZoek2IZtLjSxYhJwF2tiEPSRpIW\npP3dEhGLWhDHt4FPA2uKlrXiMwngZkn3SDqxBXGUuxtg8ybHUGoqWc2RZsYREc8D3wIeJau1roiI\nm5oZQztodqLOzS0mkf2brRZPw2KVNJzs69EZsf79l02LIyLWRMTeZLXav5V0cDPjkPQ+4OmImA+U\nvY++iZ/JARExCXg3cIqyBxqaGcfGwD7A9yJiH+AV0kWmJsawlqRNgcOBn2ywk8E/LnYG/g/Zt98x\nwHBJxzYzhnbQ7ET9OOvfLziW9f87DrblkrYFkLQd8HSFuHZIywZM0iZkSfqSiLi2VXEURMRK4Drg\nb5ocx/7AEZIeJqu5vVPSJU2OAYCIeDL9fAa4hqwjsWbGUe5ugH2Ap1p0XLwb+F36e0Bz/xZvBe6I\niOciewbjZ2RNpK36W+RSsxP12qcW03/xo8nuIWyWWWQXKkg/ry1aPlXSppLGA7sA8wa6M0kCLgQW\nR8S5LYxjdOGquaTNyNoA5zczjoj4bESMjYjxZF+zfx0RH21mDACShin1q5CaGw4luyjUzL9FpbsB\nZjcrhhLHsK7Zo7C/ZsXRC+wnabN0vhxC9qBcq/4W+dTsRnGy/973k12tnT6I+ym9Uns82YWLm4EH\ngBuBUUXlP5ti6gUOa1AMB5K1xy4gS4zzyfrpbnYcbwZ+n+L4A/DptLypcRS99ztYd9dHs/8W41n3\njMN9hWOwBXFscDdAKz4PYHOyx59HFC1r9t/iM2T/qBaS3Q2zSaNjoM0vJvrJRDPreKryZGKpg8lf\nf0QeM9HMLOecqM3Mcs6J2sws55yozcxyzonazCznnKjNzHLOidrMLOecqM3Mcs6J2sws55yozcyq\nUPnRomZIWlY0Os6UKtuvN6pRWnZUGtXmdUn79BWDE7WZWXUXk/XRUyyAf4+ISel1fZXtzyDraKq4\nv46FwAeA22oJwInazKyKiJgLvFBmVZ/9gUjaAXgP8MPi8hHRGxEP1BqDE7WZWf+cJuleSRcWD75b\notyoRnVzojYzq9/3ybrM3ZtsUOBvlRaoZVSjWm08kI3NzNrFwX0XqVlEFEacQdIPyQY6KFUY1eg9\nwFBgC0k/ioiP1bs/16jNrONFhOp59fV+aXiwgg+QXRws3We5UY3KJek+9+dEbWZWhaTLgTuA3SQ9\nJukE4OuS/iDpXrIRi85MZcdIuq7CW62960PSByQ9BuwHXCfpl1Vj8AgvZmb55hq1mVnOOVGbmeWc\nE7WZWc45UZuZ5ZwTtZlZzjlRm5nlnBO1mVnOOVGbmeXc/wdkBPcMaO7joQAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x1141dfc90>"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "###Let's make an animation!"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# for the class, we need to determine where ffmpeg lives\n",
      "# in general, you probably won't need to bother\n",
      "import platform\n",
      "from subprocess import check_output\n",
      "if 'Windows' in platform.platform():\n",
      "    if '64bit' in platform.architecture()[0]:\n",
      "        ffmpeg_path = os.path.join(binpath, 'ffmpeg.exe')\n",
      "    else:\n",
      "        ffmpeg_path = os.path.join(binpath, 'win32', 'ffmpeg.exe')\n",
      "else:\n",
      "    #Assume it is in path on macos\n",
      "    ffmpeg_path = 'ffmpeg'\n",
      "print 'ffmpeg_path is: ', ffmpeg_path\n",
      "\n",
      "figfiles = ['s_only%d.png','s%d.png','spt%d.png']\n",
      "anmfiles = ['Theis_movie1.mp4','Theis_movie2.mp4','Theis_movie3.mp4']\n",
      "\n",
      "# note the tricky way we can iterate over the elements of \n",
      "# two lists in pairs using zip\n",
      "for figfile,anmfile in zip(figfiles,anmfiles):\n",
      "    try:\n",
      "        os.remove(anmfile)\n",
      "        print 'Deleted the existing animation: ', anmfile\n",
      "    except:\n",
      "        pass\n",
      "    # now we do a system call, making the movie using command line arguments\n",
      "    # for ffmpeg\n",
      "    output = check_output([ffmpeg_path, \n",
      "                '-f', 'image2', \n",
      "                '-i', figfile, \n",
      "                '-vcodec', 'libx264',\n",
      "                '-pix_fmt', 'yuv420p',\n",
      "                anmfile])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "ffmpeg_path is:  ffmpeg\n",
        "Deleted the existing animation:  Theis_movie1.mp4\n",
        "Deleted the existing animation: "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " Theis_movie2.mp4\n",
        "Deleted the existing animation: "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " Theis_movie3.mp4\n"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# finally, we can embed the movie in the notebook\n",
      "# we could also look at it outside of the notebook with no trouble\n",
      "from IPython.display import HTML\n",
      "from base64 import b64encode\n",
      "video = open(anmfiles[2], \"rb\").read()\n",
      "video_encoded = b64encode(video)\n",
      "video_tag = '<video controls alt=\"test\" src=\"data:video/x-m4v;base64,{0}\">'.format(video_encoded)\n",
      "HTML(data=video_tag)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "<video controls alt=\"test\" src=\"data:video/x-m4v;base64,AAAAIGZ0eXBpc29tAAACAGlzb21pc28yYXZjMW1wNDEAAAAIZnJlZQAAEm9tZGF0AAACoQYF//+d3EXpvebZSLeWLNgg2SPu73gyNjQgLSBjb3JlIDE0MiAtIEguMjY0L01QRUctNCBBVkMgY29kZWMgLSBDb3B5bGVmdCAyMDAzLTIwMTQgLSBodHRwOi8vd3d3LnZpZGVvbGFuLm9yZy94MjY0Lmh0bWwgLSBvcHRpb25zOiBjYWJhYz0xIHJlZj0zIGRlYmxvY2s9MTowOjAgYW5hbHlzZT0weDM6MHgxMTMgbWU9aGV4IHN1Ym1lPTcgcHN5PTEgcHN5X3JkPTEuMDA6MC4wMCBtaXhlZF9yZWY9MSBtZV9yYW5nZT0xNiBjaHJvbWFfbWU9MSB0cmVsbGlzPTEgOHg4ZGN0PTEgY3FtPTAgZGVhZHpvbmU9MjEsMTEgZmFzdF9wc2tpcD0xIGNocm9tYV9xcF9vZmZzZXQ9LTIgdGhyZWFkcz0xMiBsb29rYWhlYWRfdGhyZWFkcz0yIHNsaWNlZF90aHJlYWRzPTAgbnI9MCBkZWNpbWF0ZT0xIGludGVybGFjZWQ9MCBibHVyYXlfY29tcGF0PTAgY29uc3RyYWluZWRfaW50cmE9MCBiZnJhbWVzPTMgYl9weXJhbWlkPTIgYl9hZGFwdD0xIGJfYmlhcz0wIGRpcmVjdD0xIHdlaWdodGI9MSBvcGVuX2dvcD0wIHdlaWdodHA9MiBrZXlpbnQ9MjUwIGtleWludF9taW49MjUgc2NlbmVjdXQ9NDAgaW50cmFfcmVmcmVzaD0wIHJjX2xvb2thaGVhZD00MCByYz1jcmYgbWJ0cmVlPTEgY3JmPTIzLjAgcWNvbXA9MC42MCBxcG1pbj0wIHFwbWF4PTY5IHFwc3RlcD00IGlwX3JhdGlvPTEuNDAgYXE9MToxLjAwAIAAAAZoZYiEADf//vbw/gU2O5jQlxHN6J0zH78VuLo0N73OBq/w7v3GYBhj0hB2oFAeZ7BwnRG2XjBsZUfDIwi2kAZYHtcAtLIFq3L4aUm0zuxnDsNbV3EwfDItNGjPLTIa0xwD9hsBziuzsBnAanyCaY9L3clLZGgu2jA5qoW6dJiw/4RLYEf94iOR2jk20jwQe6AAVqsgQzSlW7AA6kR1+gZUH7PUdlEaGHIHCzUb5xqWt4W9ti8EUYbA3LywxWAAtKmAACHx2H0cQUkgbwu3PIZShYzugfT8quKTDKm+v8LgHIgXwZgXmK4dm9nOtMw8APlgkhKw/xei5KmGSQKY/O1wrCwhnqBjTQfAyJMXWRObFlhi/Ezqrh94ChHaPWO/YlfjRpFCVsOXiqw679N6hGV5r1c/LOUAjNKemZRD7NQimW1hDI+a2VEcPgmBPGe4L/cMfZ1O1y8GIDENCcBYZoxY2Rq1Wn6Jst1B+0JSoo6JVx72VSX5WgS9ME3rsV99AML1TCFOZYjlI8g+DWfeGGOo0m4GPx1oxSCam+yOmToNtrSRdFdqqEiteV+oJQGPYJGeZIJSHQVwVDhiRW8tp8oBzkUv2a8kV4Q2dm94wvmlWMnmGBQu3tqSaUv2Pn2ZZsv+qE33HqnzESuTrxfhRnjidEnIB5ymHF1VcHrqcGx3m85esStYY6Oorgf0pZpvb/TkVm4kRCw5t1sX9eTWX/BxwOekFmFvT2+rx2gWgxiBYSSFlxIgNDCtUXUoI95l7yYs+w58CrhCWNDnN1WzjNuXEIoGbQUpV7VoE6fEVRDjIsPEEXSOaN7E14zBsOE7OUyvP78jC26Qj4ICuZ+2EEI1dfPO8VIRBH2K6xN7jBAoYsZ84rHB/5X43IDyzS9n9ZIvTBeR9tgr7nkKUBguL40vW1PaT7NDqskRXX8bRmr2K1f7YCUOVuLPdJ0WWzINj1vNvj9E32IQBPFCHn3z692mXaDmVtiDLEBlYD2T4Zog5H8B77VQAACmn3ddlKQ1+qRjuapUM/qWRdDBtBSFJqyGIUbgV2gyb1FSSLabMfsQ8gGUUnJOFEuqN8vGax4r6x5BxqxYrzxtKRbWKwg0HF4Q75iXfluBblw6ux4Jh1hl5wHz4QXWTnldy8HQqXLQ6Rf4Aq6mb6NIEI/XctFmQ8CeGVqqLMkNa/o7fGqI8OjdNS70tJ6TRJn2zK4s6tNmV1xxYW/8I96tXXMZkzRHaC7TXgRKGDB4ndBlXQ5RqQmKjZfpAEfiIxOsZAqu3PEGmjRWZ/gHP7SLf1bVCsuLE6eSMaDmeoeZYikrSVUj8ZYFu7BtcezXN+/k+gj0BcqWXTk9MSuRwvDRa5K1X7rzbv38f06RAA4eyGADIi/6lreBAGpLlvsQSVwQJcNyUDbJgHH9vdeQ9fZELoQU23GH/tjm/U6t2Qf423sUKyyxM+Ppt0MvS24NbeiouJ8qPaf3JB7DvSCDUI6P4zrRMu/G/lc/UccSArGxgg2PhVmRsZLuHLk31xIW5MoklRBXcj8zD7500WDEKT1wm7/zLUua8+k9FG2NjbXRLFmBXYM+dwpviEU+vkjC6tqzqXGvpaBrmMZY4RyuawA7PYsis6dJpGeVBI/zsdL4edPIZyLeYDBEffNb56+bWlALu7ENmT7SDYF57O/7MKASzatTjSricldOLINtF8k++KHPTBL7vYWyg4T3sYb9q5vqT/XKPMQO9icH+X0laXK3pTC95HkYHkJ5ZI7mF3maDhO35VmCzv8DN4Jmd7ZRKc3jwLxo7p2YOX/rsPRu7KZdwPOQF6MAAEChFqmSCbqtzGATE+g8NaL0NWJ1RllWkMWVTNJ6xHLGNXMxNXzd3yrFDl9VdiE70aRrD17TuhorH0oSLa7FUyJMfP6PhAcPugsdRd4qD2QiFQ7I4YdRDaIcitculT4iNUUJ9BKuF45WR9JkhpSgYfTlQtpPgyb21ite/GGs87oJ5dRlXFqMktto3vkjewYFKqLLsMLQEIW5VwyOya8in3O/S6k451tO7mDDgC5Cgz7r8ACQegMg92kewLmbzocugFbqDuGoCkjSqH2J2lVZVET/Z7vBM8TH9TL4xuaNgKYZaHvciavXt+Phbs5YE14mmwlW+qzdpmp/2lhZxWfOqhZY26A3R2kp11nIO6awi0iSFpFTnBf6weHGjIEAAAH2QZokbEM//p4QElMuQH+fcLSSec94/iUKulOSPv7en77OukhexUJ/j6dBPwlX/V4bfgmfDorEF/j9TCQ3+7n9Zp/4tzz6MkGIXzrzfbpycVAACYvSHuQbZxo1IXv5a7OJGXGtyfi4dIfGAOkNA5bk0OB7RilTsF1QWiLua9fpTPPTfXCMJWha3zrmENfJIpYPgcquNYTIasxVk3MiGysz4UcQe7fYAspz21ASGIycsjFB5HPRH4Ol1J7iePe+ixEX7Dle0IpGUNapiMnidqIgBVrF3nRY5guDRjoFZaR4xlCaoPVHLUQYySU9nvIPoVlUR5LopfxI02l/hpuoGn+VKcJDTw3kIJ2xjKg0PYu+cnhvtNk+fAB5qNJdUaUqsQM/p17Qs95f7FZfh09x077Lm1rh9zySt3pPCPa6wkA7gQf6ARhJfrIpUCi/CzQcyAtxZoAxsR3KjINMC02+U9UhA3P/WDtqN4/JRse3wZr7j4c8/2nU7NPUgRrz3xXNKAd9K1B53V2vYQ0UMmZLxKN0fdgK4wksegsMhzl6zleajwIV4VGY7yPlAZ72wNKrR0LYmgaL67fcd0LvExW4jhjExKsuUdIsI9Nt9uAZRughLa0O9yaxIWE7J+wQh6Rqd/hZDQVYuUurXnCmMtOx/r1oFUKUjrujEAAAAElBnkJ4hf8Bevjeq4C9glgA5w9b4d/KYd4UcQXXNW4an5IQjTJrk02dHvqv8/RWhvugN3QAlAboeHeSIL1OoqBx5VdhIm85wtI3AAAAdQGeYXRCvwH46GD9kiLKRkvO4BGUwAc0XRkxrpqTmAEH8YTk60mdxsv4dOl+HfcX6RzbcBI8Fal6cPfxU12oHDSkuHO8HbIXgrAW1WbuW7n4+4dhogJ6T7qr+WOzE4/RS8q+IbV/OHSlluub4F9foORXOqlQQAAAADsBnmNqQr8ADzAJLYEBDwsAITnLyF+g1hq2B8I8eRhU0aL5izLFon8pqqhvp3XgO/dkHdLtpF17d93YewAAAGFBmmVJqEFomUwIZ//+nhAARzaN2jkuNdRUyoAOGksxtXCl6pn18gTDJdWnRnIzAEirn9VcHQfM5TtGCG6veqBaogeYDEx91jVPwgi8cukptpwy8IpNe08HsFs7K1WeW6p3AAAAT0GahknhClJlMCG//qeEACTf3SNqreuN+u16t+2e8+NuUyyhZtSAbUC+BeqbuZLdMmDJcYpXA7NhjxLT0FPDEcA9PD6jN1PQF/uXggb9bSEAAABvQZqnSeEOiZTAhv/+p4QAJNyZrDD746wDX6f1vPQALLq/1o51zjIqbMRFgpNdX+V76P2hmty+QjryegaI2C9254ylz1UuzXD9Iy9upjGZDC0+CuyErk/P35EQZ+IOZ3wN0aZ8VeztqsXi8E+3A+dBAAAAJEGayEnhDyZTAh3//qmWABICe8KfnXcVeABcr5uLYFme5ZL58AAAAKJBmuxJ4Q8mUwId//6plgASF3GmaidJgEHSYWSEyRTmr0fAqfAS+Z3eloHIA2CMPuK1/vssxtgd0dGnEw8zss8NMY3qdivgyJOSH3i0n4WwX9LW8yZEdDYf9dOEkA8kJ4RTSUM5ngrWGuHiD6RetVYXPMcJQnJNZmJXImtzRjybza8Os+KQi9lcXejayrjy8op2dCr+7gb83YrIlmIWEIVgQcAAAAA1QZ8KRRE8L/8AFZHpM6xwhAAftVcBQ74k3Mp9BkYBKh+iR1Kvup/ZB67uwT/SDxr+X2dn0oEAAAAbAZ8pdEK/ABz6LdWsgt2dato1K8Z6JO+lh+b8AAAAGQGfK2pCvwAdACdXfHjTyW4cAGjsC/8xTsAAAACTQZswSahBaJlMCG///qeEACLC+ciI+AK1cyU5cNnWouFCg2aN6rtVIwfrK2A3R38AgK8sHpwDImmhhngRsc2RAnuBN6MeD0vXL/HVatLecPMs2eUym8WW4cjzg1/Ryk64o8WgNDtuOll5MgTuKTFzkMez21iaOec/1mPSbJW7/qkZ7pU1PyAsgVibU6XCvMobyBBxAAAAKUGfTkURLC//ABT8ynTgssAHZfI23MYmja+STfzPmTtJ/Ft9Y2HZlR9TAAAAHQGfbXRCvwAb8ckyIFBjYAM77kj74NHarKW7Jc+ZAAAAFAGfb2pCvwAcUCJSMfnFdaDw8xTsAAAAoEGbdEmoQWyZTAhv//6nhAAhnXSOuogBYo1GnpADoY1fhGVSw588Xhq2PQl5Uj5VkgK4xOyRJyFhTnFjlRCgln+zuKKj1IjRPmjNkTBZOH8c+k3NcEv3qsXiF9XZZ5mno6aEhjOKQfACaDq6habR5spDswAgoEQX6GX/amngONxWJI0pV1vgom0e9CwVXmdn0vXzeNhXoAzsQpYKRatbNAgAAAAsQZ+SRRUsL/8AJ9POqTZI+LBBHKi6QAB/NJWmChvYiC+GLlJ/GzJnZj2AFbEAAAAlAZ+xdEK/ADYR4wzhRVAAQ/Oa6wffCPxXN/Yym4Hfg1s8aqtLuAAAABcBn7NqQr8ANhYaYocKq5VBwoYAeH5vwAAAACVBm7VJqEFsmUwIb//+p4QAILJY97vy/IAHWC17UXDtvF8qEcHHAAAANEGb2UnhClJlMCG//qeEAD9i7OpAQOs4DioZw6AIlnJ0+u0cpPzzcHW/Wj50qBsXfeS1B6wAAAAnQZ/3RTRML/8AJ1Sds2V2HrHyRTYABLhkAK83LNCa+bi474m6YGVBAAAAHwGeFnRCvwA1/88JdDz1JXdAAh0icdj0lLeQ7SeD7ZkAAAAlAZ4YakK/ADYOmV/EWW1dtqACGEVMLoETWyriZu6x/EkJ/X062AAAAFRBmh1JqEFomUwIX//+jLAA9AWwAAugxVt4RKZR6mDYJlQVJ2UnXDIByF3IyE364NlfKs7cGPfcyGDO4kdJyRAYcbiP6LhtTYeFMmm4t4IsXfS0ZOUAAAAsQZ47RREsL/8AJ9kUjp4B4wNJNphFmkrjJgeRQAT0YLqaQ1Nv8+jaPuZIAZUAAAAkAZ5adEK/ADX/8c32JdHg80AAh/VgrljmYx5H9OK0NQ39K1NxAAAAIQGeXGpCvwA2DoZFvDrVTcvDABAR2ITVqUNQjU6i+kResQAAACNBml5JqEFsmUwIV//+OEADt8GJUantuqN+6ABAD89yfjkHLAAABG9tb292AAAAbG12aGQAAAAAAAAAAAAAAAAAAAPoAAAE2AABAAABAAAAAAAAAAAAAAAAAQAAAAAAAAAAAAAAAAAAAAEAAAAAAAAAAAAAAAAAAEAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACAAADmXRyYWsAAABcdGtoZAAAAAMAAAAAAAAAAAAAAAEAAAAAAAAE2AAAAAAAAAAAAAAAAAAAAAAAAQAAAAAAAAAAAAAAAAAAAAEAAAAAAAAAAAAAAAAAAEAAAAABsAAAASAAAAAAACRlZHRzAAAAHGVsc3QAAAAAAAAAAQAABNgAAAQAAAEAAAAAAxFtZGlhAAAAIG1kaGQAAAAAAAAAAAAAAAAAADIAAAA+AFXEAAAAAAAtaGRscgAAAAAAAAAAdmlkZQAAAAAAAAAAAAAAAFZpZGVvSGFuZGxlcgAAAAK8bWluZgAAABR2bWhkAAAAAQAAAAAAAAAAAAAAJGRpbmYAAAAcZHJlZgAAAAAAAAABAAAADHVybCAAAAABAAACfHN0YmwAAACYc3RzZAAAAAAAAAABAAAAiGF2YzEAAAAAAAAAAQAAAAAAAAAAAAAAAAAAAAABsAEgAEgAAABIAAAAAAAAAAEAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAY//8AAAAyYXZjQwFkABX/4QAZZ2QAFazZQbCWwEQAAAMABAAAAwDIPFi2WAEABmjr48siwAAAABhzdHRzAAAAAAAAAAEAAAAfAAACAAAAABRzdHNzAAAAAAAAAAEAAAABAAAA8GN0dHMAAAAAAAAAHAAAAAEAAAQAAAAAAQAACgAAAAABAAAEAAAAAAEAAAAAAAAAAQAAAgAAAAAEAAAEAAAAAAEAAAoAAAAAAQAABAAAAAABAAAAAAAAAAEAAAIAAAAAAQAACgAAAAABAAAEAAAAAAEAAAAAAAAAAQAAAgAAAAABAAAKAAAAAAEAAAQAAAAAAQAAAAAAAAABAAACAAAAAAEAAAQAAAAAAQAACgAAAAABAAAEAAAAAAEAAAAAAAAAAQAAAgAAAAABAAAKAAAAAAEAAAQAAAAAAQAAAAAAAAABAAACAAAAAAEAAAQAAAAAHHN0c2MAAAAAAAAAAQAAAAEAAAAfAAAAAQAAAJBzdHN6AAAAAAAAAAAAAAAfAAAJEQAAAfoAAABNAAAAeQAAAD8AAABlAAAAUwAAAHMAAAAoAAAApgAAADkAAAAfAAAAHQAAAJcAAAAtAAAAIQAAABgAAACkAAAAMAAAACkAAAAbAAAAKQAAADgAAAArAAAAIwAAACkAAABYAAAAMAAAACgAAAAlAAAAJwAAABRzdGNvAAAAAAAAAAEAAAAwAAAAYnVkdGEAAABabWV0YQAAAAAAAAAhaGRscgAAAAAAAAAAbWRpcmFwcGwAAAAAAAAAAAAAAAAtaWxzdAAAACWpdG9vAAAAHWRhdGEAAAABAAAAAExhdmY1Ni4xNS4xMDI=\">"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 15,
       "text": [
        "<IPython.core.display.HTML at 0x11af5e290>"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    }
   ],
   "metadata": {}
  }
 ]
}