shock_tube.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# Shock problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.optimize\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['figure.figsize'] = 11, 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "class ShockTube(object):\n",
    "    @staticmethod\n",
    "    def shock_tube_function(p4, p1, p5, rho1, rho5, gamma):\n",
    "        z = (p4 / p5 - 1.)\n",
    "        c1 = np.sqrt(gamma * p1 / rho1)\n",
    "        c5 = np.sqrt(gamma * p5 / rho5)\n",
    "        gm1 = gamma - 1.\n",
    "        gp1 = gamma + 1.\n",
    "        g2 = 2. * gamma\n",
    "        fact = gm1 / g2 * (c5 / c1) * z / np.sqrt(1. + gp1 / g2 * z)\n",
    "        fact = (1. - fact) ** (g2 / gm1)\n",
    "        return p1 * fact - p4\n",
    "\n",
    "    def __init__(self, p_1, p_5, rho_1, rho_5, gamma=1.4):\n",
    "        self.gamma = gamma\n",
    "        self.p_1 = p_1\n",
    "        self.p_5 = p_5\n",
    "        self.rho_1 = rho_1\n",
    "        self.rho_5 = rho_5\n",
    "        self.u_1 = 0\n",
    "        self.u_5 = 0\n",
    "        self.p_4 = scipy.optimize.fsolve(self.shock_tube_function, self.p_1, \n",
    "                   (self.p_1, self.p_5, self.rho_1, self.rho_5, self.gamma))[0]\n",
    "        z = (self.p_4 / self.p_5 - 1.)\n",
    "        c5 = np.sqrt(gamma * self.p_5 / self.rho_5)\n",
    "\n",
    "        gm1 = self.gamma - 1.\n",
    "        gp1 = self.gamma + 1.\n",
    "        gmfac1 = 0.5 * gm1 / self.gamma\n",
    "        gmfac2 = 0.5 * gp1 / self.gamma\n",
    "        fact = np.sqrt(1. + gmfac2 * z)\n",
    "        self.u_4 = c5 * z / (self.gamma * fact)\n",
    "        self.rho_4 = self.rho_5 * (1. + gmfac2 * z) / (1. + gmfac1 * z)\n",
    "\n",
    "        # shock speed\n",
    "        self.w = c5 * fact\n",
    "        \n",
    "        self.p_3 = self.p_4\n",
    "        self.u_3 = self.u_4\n",
    "        self.rho_3 = self.rho_1 * (self.p_3 / self.p_1)**(1. / self.gamma)\n",
    "    \n",
    "    def get_positions(self, t):\n",
    "        c1 = np.sqrt(self.gamma * self.p_1 / self.rho_1)\n",
    "        c3 = np.sqrt(self.gamma * self.p_3 / self.rho_3)\n",
    "        return np.array([- c1 * t, (self.u_3 - c3) * t, self.u_3 * t, self.w * t ,])\n",
    "    \n",
    "    def get_points(self, t, num_points=10):\n",
    "        jump = 10e-10\n",
    "        c1 = np.sqrt(self.gamma * self.p_1 / self.rho_1)\n",
    "        c5 = np.sqrt(self.gamma * self.p_5 / self.rho_5)\n",
    "        gm1 = self.gamma - 1.\n",
    "        gp1 = self.gamma + 1.\n",
    "        pos = self.get_positions(t)\n",
    "        x_pos_1_2 = np.linspace(pos[0], pos[1], num_points)\n",
    "        x_min = pos[0] - (pos[-1] - pos[0]) * 0.4\n",
    "        x_max = pos[-1] + (pos[-1] - pos[0]) * 0.4\n",
    "\n",
    "        u_1_2 = 2. / gp1 * (c1 + x_pos_1_2 / t)\n",
    "        fact = 1. - 0.5 * gm1 * u_1_2 / c1\n",
    "        rho_1_2 = self.rho_1 * fact ** (2. / gm1)\n",
    "        p_1_2 = self.p_1 * fact ** (2. * self.gamma / gm1)\n",
    "\n",
    "        x_values = np.array([x_min] + list(x_pos_1_2) + \n",
    "                            [pos[2], pos[2] + jump, pos[3], pos[3] + jump, x_max])\n",
    "        p_values = np.array([self.p_1] + list(p_1_2) +\n",
    "                            [self.p_3, self.p_4, self.p_4, self.p_5, self.p_5])\n",
    "        u_values = np.array([self.u_1] + list(u_1_2) +\n",
    "                            [self.u_3, self.u_4, self.u_4, self.u_5, self.u_5])\n",
    "        rho_values = np.array([self.rho_1] + list(rho_1_2) +\n",
    "                            [self.rho_3, self.rho_4, self.rho_4, \n",
    "                             self.rho_5, self.rho_5])\n",
    "        return x_values, p_values, rho_values, u_values, x_min, x_max"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAo4AAAFACAYAAAAyBP25AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl0VdXdxvHvTkgYJSCTA0lAGZRBIiKDTLYFywxaRHAC\nAZFUE5AhoJUp2paAzJWglalaUbSWKYCClkEEFWVGsWghxIlBAREhSPb7x4lvUwzkJLk3596b57MW\nK7n37nP203Wa5Oc+Z+9trLWIiIiIiOQlzOsAIiIiIhIcVDiKiIiIiCsqHEVERETEFRWOIiIiIuKK\nCkcRERERcUWFo4iIiIi4osJRRERERFxR4SgiIiIirqhwFBERERFXSnjVceXKlW2NGjW86l5ERERE\nsn344YdHrbVV8mrnWeFYo0YNtm7d6lX3IiIiIpLNGHPQTTvdqhYRERERV1Q4ioiIiIgrKhxFRERE\nxBXPnnEUERERCQXnzp0jIyODM2fOeB0lT6VKlaJ69epEREQU6HgVjiIiIiKFkJGRwWWXXUaNGjUw\nxngd56KstRw7doyMjAxq1qxZoHPoVrWIiIhIIZw5c4ZKlSoFdNEIYIyhUqVKhRoZVeEoIiIiUkiB\nXjT+rLA5VTiKiIiIiCsqHEVERETEFRWOIiIiIuKKZlWLiIiIhIDevXtjreXAgQN8/fXXzJ49m86d\nO/u0DxWOIiIiIj4ydChs3+7bc8bFwfTpebfbsWMHPXr04JVXXuGdd95h2LBhPi8cdataREREJMj9\n+OOPHD16lHHjxgFQr149vvvuO5/3oxFHERERER9xMzLoD7t376Z27dqUKlUKgI8++ohGjRr5vB8V\njiIiIiJBbseOHaSnp3PmzBnOnz/PuHHjmDRpks/7UeEoIiIiEuR27NjBPffcw6233srJkyd5/PHH\nadmypc/7UeEoIiIiEuR27NjBX//6V1JSUvzaT56TY4wx84wxh40xuy/yuTHGzDTG7DfG7DTGNPZ9\nTBERERG5mM8++4zatWv7vR83s6oXAB0u8XlHoHb2v0FAauFjiYiIiIhbX3zxBWFh/l8sJ89b1dba\nDcaYGpdo0h34m7XWAluMMRWMMVdaa7+61HnPnT7H1x99ma+wUnyUrlyWqJgor2OIiIhIDr54xvFq\n4FCO1xnZ712ycIz4eCdX3HS1D7qXUJRJBHv+9Br1H+vmdRQRERHJ5ovC0eTyns21oTGDcG5nE1um\nGhtuT/ZB9xKKKi6eQ+yYvmTdtY2wa2p4HUdERETwTeGYAUTneF0dyPUetLX2OeA5gCZNmtg2Lw7y\nQfcSil6Pa0f0yMZ8d9tdVNq7ESIjvY4kIiJS7PniKcplwP3Zs6ubAyfyer5RJC/dH72GCTXmU+mz\n9/lpeJLXcURERAR3y/EsAjYDdY0xGcaYAcaYwcaYwdlNVgKfA/uBvwK/91taKTbCw6Hjs7cznSGU\n+MsMeP11ryOJiIgUe25mVffJ43MLPOyzRCLZbrsNZrSbROt/babxA/0xcXFwzTVexxIRESm2/L/g\nj0gh/HFyJD3Pv8KPZw306gVnz3odSUREpNhS4SgBLS4O2txfg/vOL4QPP4Thw72OJCIiErBmz55N\ngwYNiI2NZdasWT4/v/aqloD35JNQ55VurKoznI7PTIE2bZzRRxERkUAzdChs3+7bc8bFwfTpeTb7\nxz/+wZo1a9i2bRtHjx6lYcOGxMfHU6KE78o9jThKwIuJcX4Ou+/9Mz/c0AIGDoR//9vrWCIiIgFl\n5syZpKSkEBERwZVXXklERARZWVk+7UMjjhIURo+Gv/41gkFRr/D3iDhnxHHzZihVyutoIiIi/+Vi\nZNAfzp07x86dO6lTpw4AX331FZUrVybSx+sga8RRgkKFCjBmDLy0MZqPHn3BuQ0wdKjXsURERALC\n3r17OXHiBJ9//jlZWVk89thjJCQk+LwfFY4SNH7/e2c1nv6vdSJr5Ch49llYtMjrWCIiIp7btm0b\n99xzD3369OGGG24gJiaGQYN8v0OfblVL0IiMhD/9CXr3hhcTn+L+Vptg0CBo3Bjq1vU6noiIiGe2\nb99Oly5duOuuu/zaj0YcJajceSfcfDP8YVwJfpz/svOM4513wunTXkcTERHxzPbt24mLi/N7Pyoc\nJaiEhcHkyZCRATP/cTW8+CLs3g2JiV5HExER8cy6deuoWwR333SrWoJO27bQtatz23rAZ7+l8uOP\nwx//6Hxw331exxMfyLJZJK9PZs+RPV5HEQkJ5SPLM6PjDMpFlvM6igQ5FY4SlCZOhIYN4amnYPrT\n4+Gdd2DwYLjpJqhXz+t4UkgLty9kwvoJ1Lq8FpHhvl1KQqS4OZV5ivQT6QxoPIBbom/xOo4EORWO\nEpTq1YMBA2D2bEhIKMG1ixY5K+vfeSe8/z6ULet1RCmgY6ePMXLNSFpGt2TDAxsIM3qiRqQw3vzs\nTX774m+9jiEhQr+RJWhNmAAREfCHPwBXXgkvvQQff+zMtLbW63hSQKPWjuL4meOkdk5V0SgiEmD0\nW1mC1pVXwogR8MorziAjv/mNs7H1Sy/B1Klex5MC2JS+ibnb5jKsxTAaVmvodRwREbmACkcJaiNG\nQNWqMHJk9iDj449Dz56QlARvvOF1PMmHc+fPMThtMNHloxnbdqzXcUREJBcqHCWoXXYZjB8PGzbA\nihWAMTB/PtSv76wUvn+/1xHFpRnvzWD34d3M6jhLMz9FRAKUCkcJegMHQp06MGoU/PQTUK4cLF3q\nLPrYvTt8/73XESUP6SfSGbduHF3rdKX7dd29jiMiIhehwlGCXkQEpKQ482Lmzct+s2ZNePVV2LcP\n7r8fsrI8zSiXNmT1EABmdZzlcRIRkeDVvHlzDhw4AMAXX3xBkyZNfN6HluORkNC9O7RsCePGwd13\nO4OO/PrXMGUKDB0KycnOPW0JOMv3LWfJJ0uY+JuJxFaI9TqOiEihDF09lO1fb/fpOeOuiGN6h+mX\nbGOtJT09ndhY5/fozp07adjQ95MMNeIoIcEYZyvCr792vv6/xETo29dZu+ef//Qsn+Tuh8wfSFiV\nQP0q9RnWYpjXcUREgtb+/fupWbMmxhjAf4WjRhwlZLRo4cyHSUlx7k5fey1ORTlnjnMf+/77YcsW\nZ+KMBIQnNzzJwRMH2dBvAxHhEV7HEREptLxGBv1l165d/1Mobt26lYceesjn/WjEUULKlCkQGQmP\nPJJjDfBSpZzRxnLlnHva337raUZx7Dm8hymbp/BA3AO0jm3tdRwRkaD27bffUrp0aQA+/vhj0tLS\ndKtaJC9XXeWsAb569QV3pq+6Cl5/HQ4dcoYlf/rJs4ziPIsTnxZP+ZLlmdR+ktdxRESC3m9/+1ve\neustevXqxauvvkqlSpWoVq2az/tR4Sgh5+GHoVEjGDIETp3K8UGLFs7m1mvWwOjRnuUTWLhjIRvT\nNzKp3SQql6nsdRwRkaAXHR3Nzp07Wbx4MWPHjuXQoUN+6UeFo4ScEiUgNRUyMpzJ1P9jwADnPvaU\nKfD3v3uSr7g7dvoYI94cQcvoljxw4wNexxERkXxQ4SghqUULZ2HwadNg9+4LPpw6FW691WmwdasX\n8Yq10WtHc/zMcVI7pxJm9CtIRCSY6Le2hKyJEyEqCn7/+xwTZcBZMXzxYqhWDW6/Hb75xrOMxc2m\n9E08v+15hrUYRsNqvn9oW0TEK/Z//tAErsLmVOEoIatSJWdpno0b4YUXLviwShVYsgSOHYPf/Q4y\nMz3JWJycO3+O+LR4ostHM7btWK/jiIj4TKlSpTh27FjAF4/WWo4dO0apUqUKfA6t4ygh7YEHYO5c\nGDECunaFihVzfBgXB/PnO7OsExOd9R7Fb2a8N4Ndh3ex5K4llIss53UcERGfqV69OhkZGRw5csTr\nKHkqVaoU1atXL/DxKhwlpIWFORNlGjeGP/zBmVT9P+66C7Zvd+5rx8XB4MGe5Ax16SfSGbduHF3r\ndKX7dd29jiMi4lMRERHUrFnT6xhFQreqJeQ1avTfAcUPPsilwVNPQadOkJAAa9cWeb7iYMjqIQDM\n6jjL4yQiIlIYKhylWJgwAa64AuLj4fz5Cz4MD4dFi+D6653nHXft8iRjqFq+bzlLPlnC2DZjia0Q\n63UcEREpBBWOUiyUL+8szfPhh/DssxdpsHKlsy1hp07w5ZdFnjEU/ZD5AwmrEqhfpT7DWgzzOo6I\niBSSCkcpNnr1gnbt4PHH4euvc2lQvTqkpcHx49C5M3z/fZFnDDVPbniSgycOkto5lYjwCK/jiIhI\nIalwlGLDGHjmGfjxRxg58iKN4uLg1Ved29W9emlP60LYc3gPUzZP4YG4B2gd29rrOCIi4gMqHKVY\nqVMHkpLgxRdh3bqLNOrQwZmKvXq1s/F1gK/LFYistcSnxVO+ZHkmtZ/kdRwREfERV4WjMaaDMWaf\nMWa/MWZ0Lp/HGGP+ZYzZZozZaYzp5PuoIr7x+ONQs6azo8xF1/1+8EF47DF47jmYpMInvxbuWMjG\n9I1MajeJymUqex1HRER8JM/C0RgTDjwDdATqAX2MMfUuaPYEsNhaeyPQG7hwtTyRgFG6NMycCR9/\n7EyYuainnoI+fWD0aHj55SLLF+yOnT7GiDdH0DK6JQ/c+IDXcURExIfcjDg2BfZbaz+31mYCLwMX\nruBrgfLZ30cBmpIqAa1LF+jRA5KT4eDBizQKC3N2lmndGvr2dfYulDyNXjua42eOk9o5lTCjp2FE\nREKJm9/qVwOHcrzOyH4vp/HAvcaYDGAlkJDbiYwxg4wxW40xW4NhWx4JbdOnO1+HDr1Eo5IlnT2t\na9Z0Ks19+4okW7DalL6J57c9z7AWw2hYraHXcURExMfcFI4ml/cunC3QB1hgra0OdAJeMOaXQw3W\n2uestU2stU2qVKmS/7QiPhQbC2PHOnXhihWXaHj55c4aj+HhzhqPhw8XWcZgcu78OQanDSYmKoZx\nbcd5HUdERPzATeGYAUTneF2dX96KHgAsBrDWbgZKAXoiXgLeo486G8YkJMDp05doeM01sHw5fPUV\ndOuWR+PiacZ7M9h9eDczO8ykbGRZr+OIiIgfuCkcPwBqG2NqGmMicSa/LLugTTrwGwBjzPU4haPu\nRUvAi4yE2bPhwAH485/zaNysGbz0Erz/Ptx7by57FxZf6SfSGbduHN3qdqP7dRc+Ai0iIqEiz8LR\nWvsT8AjwBvAxzuzpPcaYZGNMt+xmw4EHjTE7gEVAP2u1+J0Eh1tvderAlBTYuTOPxj16OFOx//nP\nS6wiXvwMWT0EgJkdZnqcRERE/KmEm0bW2pU4k15yvjc2x/d7gZa+jSZSdKZNgzVr4P77nQHFyMhL\nNB4yBP7zH+egmjWd+9zF2LJ9y1jyyRJS2qUQWyHW6zgiIuJHWitDBKhc2Vnre8cOePJJFwdMmeKM\nPg4ZAkuX+j1foPoh8wcSViVQv0p9Hm3+qNdxRETEz1Q4imTr1g369XOedXz//Twah4fD3/8ON9/s\nLBK+ZUtRRAw4T254kvQT6aR2TiUiPMLrOCIi4mcqHEVymD4drrrKWe/7xx/zaFymjDPT+qqroGNH\n2LWrSDIGit2HdzNl8xT6x/WndWxrr+OIiEgRUOEokkNUFMybB598Ak884eKAqlVh7VooWxZuuw32\n7/d7xkCQZbOIT4unfMnypLRP8TqOiIgUERWOIhdo1w5+/3tn7ourXQZr1HBm1pw75xyckeHviJ5b\nuH0h76S/w+T2k6lcRku2iogUFyocRXKRkuJMmO7XD06dcnHA9dfDG2/At99C+/YQwltqHjt9jJFr\nRtIyuiX94vp5HUdERIqQCkeRXJQrBwsWOKvuJCW5POimm5y9Cw8ccJ55PHnSjwm9M2rtKE6cPUFq\n51TCfrmzqIiIhDD91he5iNatYdgwSE2FN990eVCbNvDaa866Pl27uphhE1w2pW9i7ra5DGs+jIbV\nGnodR0REipgKR5FLeOop5y70gAFw/LjLgzp3hhdecB6Q7NkTMjP9mrGonDt/jsFpg4mJimFs27F5\nHyAiIiFHhaPIJZQqBQsXwldfwdCh+Tiwd2+YMwdWrnS2owmBfa2nb5nO7sO7mdVxFmUjy3odR0RE\nPKDCUSQPN98Mjz/uFJD52iRm0CBnls0rrzjTtIN4+/b0E+mMXz+ebnW70a1ut7wPEBGRkKTCUcSF\nJ56AuDinFjx6NB8HJiXBY485+xk+9pjf8vlb4qpEAGZ2mOlxEhER8ZIKRxEXIiPhb3+D776Dhx/O\n58F//CPExzujjxMn+iWfPy3bt4yl+5Yyru04YivEeh1HREQ8pMJRxKWGDWHCBFi82Ln77Jox8Je/\nwN13O6OOqal+y+hrP2T+QMKqBOpXqc+jzR/1Oo6IiHhMhaNIPowcCc2aOY8sfvVVPg4MC3MWhuzS\nxRmyfOklf0X0qeT1yaSfSGdOlzlEhEd4HUdERDymwlEkH0qUcCbJnD7tPO+Yr/kuERHOcGWbNs5M\n6+XL/ZbTF3Yf3s3ULVPpH9efVjGtvI4jIiIBQIWjSD7Vres8qrhihTOImC+lS8OyZXDjjXDnnfCv\nf/kjYqFl2Szi0+KJKhlFSvsUr+OIiEiAUOEoUgAJCdC2LQwZAgcP5vPg8uVh1Sq49lrn1vW6df6I\nWCgLty/knfR3mNR+EpXLVPY6joiIBAgVjiIFEBYG8+c7t6r794esrHyeoHJlePttqFEDOnVyvg8Q\nx04fY+SakbSKaUW/uH5exxERkQCiwlGkgGrWhClTnJpv2rQCnKBaNedW9TXXONsUrl3r84wFMWrt\nKE6cPUFq51TCjH5FiIjIf+mvgkghPPgg3H47jB4N775bgBNUreoUj7VrQ9eusGaNzzPmx6b0Tczd\nNpdhzYfRoGoDT7OIiEjgUeEoUgjGwLx5EBMDd92Vz11lflalCrz1FtSp4xSPb7zh85xunDt/jsFp\ng4mJimFs27GeZBARkcCmwlGkkCpUgFdfhcOH4b77CvC8I/y3eLz+euje3Zk8U8Smb5nO7sO7mdVx\nFmUjyxZ5/yIiEvhUOIr4QOPGMGMGrF5diF0FK1d2nnOsVw969ICVK32a8VIOHj/I+PXj6Va3G93q\ndiuyfkVEJLiocBTxkYcegj59YMyYQizPWKmSUzw2bOg8PLlihU8zXsyQ1UMAmNlhZpH0JyIiwUmF\no4iPGAPPPec8qtinD3z9dQFPdPnlziSZG26AO+5wFgz3o2X7lrF031LGtR1HbIVYv/YlIiLBTYWj\niA+VK+c873jypFM8nj9fwBNVrOgUj3Fx0LMnLFni05w/+yHzBxJWJVC/Sn0ebf6oX/oQEZHQocJR\nxMcaNIDUVGdDmPHjC3GiChWc4rFxY2d7wtdf91HC/0pen0z6iXTmdJlDRHiEz88vIiKhRYWjiB/0\n7evsKPPUU86EmQKLioI334Sbb4ZeveC113yWcffh3UzdMpX+cf1pFdPKZ+cVEZHQpcJRxE9mzXLm\nuNx7Lxw6VIgTlS/vVJ/NmkHv3s698ELKslnEp8UTVTKKlPYphT6fiIgUDyocRfykTBmnxjt71qn3\nzp0rxMl+Lh5btHAennzllUJlW7h9Ie+kv8Ok9pOoXKZyoc4lIiLFhwpHET+qWxeef97ZjvCxxwp5\nsssucxYGv+UWuPtu58QFcOz0MUauGUmrmFb0i+tXyFAiIlKcqHAU8bO77oKHH4YpU2Dp0kKerFw5\np3i87TZno+yUFLA2X6cYtXYUJ86eILVzKmFGvwJERMQ9/dUQKQJTpsBNNzmTZj7/vJAnK1vWqUD7\n9IHRo2HkSNf7HG5K38TcbXMZ1nwYDao2KGQQEREpblQ4ihSBkiX/O6elVy/nucdCiYyEF1+ERx5x\nqtL+/fN8iPLc+XMMThtMTFQMY9uOLWQAEREpjlQ4ihSRmjVh4UL48EMYPtwHJwwLg5kzYcIE58R3\n3AE//njR5tO3TGf34d3M7DCTspFlfRBARESKG1eFozGmgzFmnzFmvzFm9EXa9DLG7DXG7DHGvOTb\nmCKhoXt3p2h85plCT4x2GANjx8Ls2ZCW5jz7ePz4L5odPH6Q8evH061uN7pf190HHYuISHGUZ+Fo\njAkHngE6AvWAPsaYehe0qQ08BrS01tYHhvohq0hI+POfnVV1Bg6ETz/10Unj4+Hll+G996BtW/jq\nq//5eMjqIQDM7DDTRx2KiEhx5GbEsSmw31r7ubU2E3gZuHDI4kHgGWvtdwDW2sO+jSkSOiIinNHG\nkiWhR49cBwgLplcvZ9Txs8+gVSvnK7Bs3zKW7lvK+Lbjia0Q66PORESkOHJTOF4N5Nz3IiP7vZzq\nAHWMMZuMMVuMMR1yO5ExZpAxZqsxZuuRI0cKllgkBERHO7sH7t/v1HuFWhw8p/bt4e234cQJaNmS\nHz7cQsKqBBpUbcDQ5roRICIiheOmcDS5vHfhwnElgNrArUAf4HljTIVfHGTtc9baJtbaJlWqVMlv\nVpGQcuut8OyzsGYNJCbmeznGi2vaFDZuhIgIkse1Jf1EOqmdU4kIj/BRByIiUly5KRwzgOgcr6sD\nX+bSZqm19py19j/APpxCUkQu4YEHYNQomDMHZszw4Ymvv57dafOZelMm/XeE02rHdz48uYiIFFdu\nCscPgNrGmJrGmEigN7DsgjZLgF8BGGMq49y6LuwyxyLFwp/+5KykM2wYLF/um3Nm2SziP5xAVJnL\nSfmmIdx+OyxY4JuTi4hIsZVn4Wit/Ql4BHgD+BhYbK3dY4xJNsZ0y272BnDMGLMX+Bcw0lp7zF+h\nRUJJWBi88IKzs0yfPrB9e+HPuWD7At5Jf4dJt02m8uoN8KtfOcObTz9d+JOLiEixVcJNI2vtSmDl\nBe+NzfG9BYZl/xORfCpTBpYtcx5P7NoV3n8frryyYOc6evooSWuSaBXTin5x/cCEwYoVcN99zvaE\nX3zhFJDh4T793yAiIqFPO8eIBIgrr3Tqu+++g27d4PTpgp1n1JpRnDh7gtTOqYSZ7B/xkiVh0SJn\nFs706c46QKdO+S68iIgUCyocRQJIo0ZOfffhh84AYVZW/o5/J/0d5m2fx7Dmw2hQtcH/fhge7szA\n+ctfYOVKaN0aMjJ8F15EREKeCkeRANO1K0yZAq+/Do8/7v64c+fPMXjFYGKiYhjbduzFGz788H8X\nCm/a1KlSRUREXFDhKBKAhg6FwYMhJQXmz3d3zLQt09hzZA+zOs6ibGTZSzfu0AHefRciI52Rx3/+\ns/ChRUQk5KlwFAlAxsDMmc5GMIMGwbp1l25/8PhBJqyfQPe63elWt9ulG/+sQQNnb+sbboDf/Q4m\nTfLhKuQiIhKKVDiKBKiICFi8GGrXdtZ5/PTTi7dNXJ0IwMyOM/PXSbVq8K9/OfsejhoFDz4ImZmF\nSC0iIqFMhaNIAKtQwZlpHR4OXbrAt9/+ss3ST5aybN8yxrcdT0xUTP47KV0aXnoJnngC5s51bmN/\np51mRETkl1Q4igS4a66BJUvg4EFn5DHngOCpzFMkrEqgQdUGDG0+tOCdhIXBk0/C3/4GmzZBixaw\nf3/hw4uISEhR4SgSBFq2dCbJrF8PDz3030cRk9cnc+jkIVI7pxIRHlH4ju67D9auhaNHoVkz2Lix\n8OcUEZGQocJRJEjcfTeMG+dsOZ2SAru+2cW0LdMYcOMAWsW08l1HrVvDli1QpQr85jfOKKSIiAgu\ntxwUkcAwbpwzSeaxx7OYSzxRJaNIaZfi+45q1YLNm6FnT+jb1+k0Odm5pS0iIsWW/gqIBBFjnFvW\n9e9ZwP6zm7izwmQqlankn84qVoTVq2HgQPjjH+Guu+D77/3Tl4iIBAUVjiJB5vvzR/mqwUjKH2/F\nXx/uy4oVfuwsIgKeew6eftrZyqZZM9i3z48diohIIFPhKBJkRq0ZxcnMk7yRkMqNcWH07AlvveXH\nDo2B4cNhzRo4cgRuvlk7zYiIFFMqHEWCyMaDG5m3fR7Dmg+j+TUNWL3aWSC8e3dnB0G/+vWv4aOP\n4PrrnXWBHnsMzp/3c6ciIhJIVDiKBIlz588RnxZPTFQMY9uOBaBSJWcg8KqroFMnp67zq+ho2LDB\n2Qdx4kRnsfCjR/3cqYiIBAoVjiJBYtqWaew5sodZHWdRNrLs/79/xRXOreoKFeC222DvXj8HKVkS\nnn3W2WVm40a46Sb44AM/dyoiIoFAhaNIEDh4/CAT1k+ge93udKvb7RefR0c763ZHREC7dvDZZ0UQ\nqn9/Z5cZY6BVK3j++SLoVEREvKTCUSQIJK5OBGBmx5kXbVOrllM8ZmY663YfOlQEwW66CbZuhbZt\n4cEHnX9nzhRBxyIi4gUVjiIBbuknS1m2bxnj244nJirmkm3r14c334TvvnNGHr/5pggCVq4Mq1bB\n4487o46tW0N6ehF0LCIiRU2Fo0gAO5V5ioRVCTSo2oChzYe6OqZxY1i5EjIyoH17+PZbP4cECA93\nFglfssTZZaZxY2f4U0REQooKR5EAlrw+mUMnDzGn8xwiwiNcH9eyJSxd6qzV3bEjnDzpx5A5de/u\nTJSpVg1++1tn5rW1RdS5iIj4mwpHkQC165tdTN08lQE3DqBlTMt8H9+uHbz2mrNET9eucPq0H0Lm\npk4deO89uPNOZ63H3/0Ojh8vos5FRMSfVDiKBKAsm0V8WjwVSlUgpV1Kgc/TtSu8+KKzas4dd8DZ\nsz4MeSnlysGiRTB1KixfDnFxzgxsEREJaiocRQLQ/G3z2XRoE5PbT6ZSmUqFOtdddzlzVt54A/r0\ngZ9+8lHIvBgDjz7qFIzh4dCmDSQna7cZEZEgpsJRJMAcPX2UpLVJtI5pTd+4vj45Z//+MGOGs8V0\nr15FOPII0LQpbNvmVK3jxjlbFxbJWkEiIuJrKhxFAkzSmiROnj1JaudUwozvfkQTE/9bPHbtCqdO\n+ezUeStf3rln/re/OQ9dNmoEr79ehAFERMQXVDiKBJCNBzcyf/t8hrcYTv2q9X1+/sREWLDA2aKw\nfXtnvccidd99zuhjrVrOpJmHHirCWTsiIlJYKhxFAsS58+eIT4snNiqWMW3G+K2fvn3h1Vedgb9b\nb4Wvv/aaI3T2AAAWvElEQVRbV7mrVQveeQeSkuC556BJE9i5s4hDiIhIQahwFAkQ07ZMY8+RPczq\nOIuykWX92tcdd0BamrOndevWcOCAX7v7pchISEmBNWucYc+mTWHWLK35KCIS4FQ4igSAA8cPMH7d\neHpc14OudbsWSZ/t2jmbuxw9Cq1awccfF0m3vwyxc6fzNTERunWDI0c8CCIiIm6ocBQJAImrEjHG\nMKPDjCLtt3lzWL/eWaKndWv48MMi7d5RpYqz1uOMGc5G240aOQ9hiohIwCnhdQCR4m7pJ0tZ/uly\nJrefTExUTJH3f8MNziOH7drBr34FK1Y4Sy4WKWOcEce2baF3b2fmTlIS3H+/85kEv4gIuPZaXU+R\nIKfCUcRDpzJPkbAqgQZVGzCk2RDPctSq5azT3b69s8X0a69B584eBGnUyBn2fPRR5xnIlILvmiMB\naNEi5z8MRCRoqXAU8VDy+mQOnTzEot8tIiI8wtMsV18NGzZAhw7Qo4ez5GKfPh4EKVMGnn3Wmf6t\nhcJDw8mTMGgQHDvmdRIRKSQVjiIe2fXNLqZunsrAGwfSMqal13EAqFwZ3n7bmaNyzz1w4gQMHuxR\nmFtu8ahj8bkjR5zCUUSCnqvJMcaYDsaYfcaY/caY0Zdo19MYY40xTXwXUST0ZNks4tPiqVi6IhPb\nTfQ6zv8oXx5WrYJOnSA+HiYGVjwREfFQnoWjMSYceAboCNQD+hhj6uXS7jIgEXjP1yFFQs38bfPZ\ndGgTk9pNolKZSl7H+YXSpZ2tCe++Gx57DEaP1hKLIiLibsSxKbDfWvu5tTYTeBnonku7J4FJwBkf\n5hMJOUdPHyVpbRKtY1rTL66f13EuKiICXnjBGXVMSXHuNJ4753UqERHxkpvC8Wog5xPqGdnv/T9j\nzI1AtLV2hQ+ziYSkpDVJnDx7ktTOqZgAX5okLAyeeQaeeAKefx5uu03zG0REijM3hWNuf9n+/6aV\nMSYMmAYMz/NExgwyxmw1xmw9ot0hpBjaeHAj87fPZ3iL4dSvWt/rOK4YA08+6Yw+bt4MzZp5tMuM\niIh4zk3hmAFE53hdHfgyx+vLgAbAOmPMAaA5sCy3CTLW2uestU2stU2qVKlS8NQiQSjzfCbxafHE\nRsUyps0Yr+Pk2733wrp1cOqUs+PMqlVeJxIRkaLmpnD8AKhtjKlpjIkEegPLfv7QWnvCWlvZWlvD\nWlsD2AJ0s9Zu9UtikSA1bfM09hzZw6yOsygbWdbrOAXSvDm8/z5ccw106QLTpmnSjIhIcZJn4Wit\n/Ql4BHgD+BhYbK3dY4xJNsZ083dAkVBw4PgBJqyfQI/retC1blev4xRKTIyzRWGPHjBsGDz4IGRm\nep1KRESKgqsFwK21K4GVF7w39iJtby18LJHQkrgqkTATxowOM7yO4hNly8Krr8L48c7zj59+Cv/4\nB+gJFBGR0OZqAXARKbilnyxl+afLGX/reGKiYryO4zNhYZCc7Gw//MEH0LQp7N7tdSoREfEnFY4i\nfnQq8xQJqxJoWLUhQ5oN8TqOX/TuDevXw9mz0KIFrNCiXCIiIUuFo4gfJa9P5tDJQ6R2TiUiPMLr\nOH7TtKkz6li3rrPP9eTJmjQjIhKKVDiK+Mmub3YxdfNUBt44kJYxLb2O43dXXw0bNsCdd0JSEjzw\ngDMKKSIioUOFo4gfZNks4tPiqVi6IhPbTfQ6TpEpUwZefhkmTICFC+FXv4JvvvE6lYiI+IoKRxE/\nmL9tPpsObWJy+8lUKlPJ6zhFyhgYO9aZdb19O9x8M2zZ4nUqERHxBRWOIj529PRRktYm0TqmNX0b\n9fU6jmd69nTWewwPh9atYdIkyMryOpWIiBSGCkcRH0tak8TJsydJ7ZyKMblt9V58NG4M27Y5i4WP\nGgUdO+rWtYhIMFPhKOJDGw9uZP72+QxvMZz6Vet7HScgVKgAixfDnDnO5Jm4OFi71utUIiJSECoc\nRXwk83wm8WnxxEbFMqbNGK/jBBRj4KGHnH2uK1aE226DP/wBfvrJ62QiIpIfKhxFfGTa5mnsObKH\nWR1nUTayrNdxAlLDhs56j/37w5/+BG3bQnq616lERMQtFY4iPnDg+AEmrJ9Aj+t60LVuV6/jBLSy\nZeH55+Gll2DXLmjUCP75T69TiYiIGyocRXwgcVUiYSaMGR1meB0laPTp40ycqVUL7rgDHnkEzpzx\nOpWIiFyKCkeRQlr6yVKWf7qc8beOJyYqxus4QeXaa2HTJhg+HJ55Bpo3h337vE4lIiIXo8JRpBBO\nZZ4iYVUCDas2ZEizIV7HCUqRkfD005CWBl98ATfd5Ow6IyIigaeE1wFEglny+mQOnTzEot8tIiI8\nwus4Qa1TJ2enmXvvhX79nCV7+vVzZmRLcIs4Dq1xRpO/eNvrNMXP9u+cr9Z6m0NCgwpHkQLa9c0u\npm6eysAbB9IypqXXcULC1Vc7BeOf/gTjx8OLL3qdSHyhMnAEmDkLZs/yOk0xdC1wH3z+ObTU0zRS\nSCocRQogy2YxOG0wFUtXZGK7iV7HCSnh4TBmjDN55ssvvU4jvhBxHOgOQ4fCXbd7nab4eeFdeP4s\nnD7tdRIJBSocRQpg/rb5vHvoXeZ3n0+lMpW8jhOSatVy/kkIOOJ8qV0LarfxNkpxtOlr4GOvU0io\n0OQYkXw6evooSWuTaB3Tmr6N+nodR0REpMiocBTJp6Q1SZw8e5LUzqkYzdwQEZFiRIWjSD5sPLiR\n+dvnM6LFCOpXre91HBERkSKlwlHEpczzmcSnxRMbFcuYtmO8jiMiIlLkNDlGxKVpm6ex58gelvdZ\nTpmIMl7HERERKXIqHEVcOHD8ABPWT6DHdT3oUqeL13FEgtO0abB4sdcpip06Ed9Cawg7/b3XUSQE\nqHAUyYO1loRVCYSZMGZ0mOF1HJHgc/nlzpZAGRleJymWws/8AECprw54mkNCgwpHkTws3beUFZ+u\nYHL7ycREadsFkXwLD4cXXvA6RbH1xR//DD897nUMCRGaHCNyCacyT5G4KpGGVRsypNkQr+OIiIh4\nSiOOIpcwYd0EDp08xKLfLSIiPMLrOCIiIp7SiKPIRez8ZifTtkxj4I0DaRnT0us4IiIinlPhKJKL\nLJtFfFo8FUtXZGK7iV7HERERCQi6VS2Si3nb5vHuoXeZ330+lcpU8jqOiIhIQNCIo8gFjvxwhFFr\nR9Emtg19G/X1Oo6IiEjAUOEocoGktUmcPHuS2Z1mY4zxOo6IiEjAUOEoksOGgxtYsH0BI1qMoH7V\n+l7HERERCSgqHEWyZZ7PJD4tntioWMa0HeN1HBERkYCjyTEi2aZunsreI3tZ3mc5ZSLKeB1HREQk\n4LgacTTGdDDG7DPG7DfGjM7l82HGmL3GmJ3GmLeMMbG+jyriPweOHyB5fTI9rutBlzpdvI4jIiIS\nkPIsHI0x4cAzQEegHtDHGFPvgmbbgCbW2huA14BJvg4q4i/WWhJWJRBmwpjRYYbXcURERAKWmxHH\npsB+a+3n1tpM4GWge84G1tp/WWtPZ7/cAlT3bUwR/1m6bykrPl3BhFsnEBMV43UcERGRgOWmcLwa\nOJTjdUb2exczAFiV2wfGmEHGmK3GmK1Hjhxxn1LET05lniJhVQINqzYksVmi13FEREQCmpvCMbeF\n7GyuDY25F2gCTM7tc2vtc9baJtbaJlWqVHGfUsRPJqybQMbJDOZ0mUNEeITXcURERAKam1nVGUB0\njtfVgS8vbGSMaQf8AWhrrT3rm3gi/rPzm51M2zKNgTcO5JboW7yOIyIiEvDcjDh+ANQ2xtQ0xkQC\nvYFlORsYY24EngW6WWsP+z6miG9l2SwGrxhMxdIVmdhuotdxREREgkKeI47W2p+MMY8AbwDhwDxr\n7R5jTDKw1Vq7DOfWdDng1ewt2tKttd38mFukUOZtm8fmjM3M7z6fSmUqeR1HREQkKLhaANxauxJY\necF7Y3N8387HuUT85sgPR0hak0Sb2Db0bdTX6zgiIiJBQ1sOSrGTtDaJ7zO/Z3an2WSPkIuIiIgL\nKhylWNlwcAMLti9gRIsR1K9a3+s4IiIiQUWFoxQbmecziU+LJzYqljFtx3gdR0REJOi4esZRJBRM\n3TyVvUf2srzPcspElPE6joiISNDRiKMUCweOHyB5fTK3X3c7Xep08TqOiIhIUFLhKCHPWkvCqgTC\nTBgzOszwOo6IiEjQ0q1qCXlL9y1lxacreLr900RHRed9gIiIiORKI44S0k5lniJhVQINqzYksVmi\n13FERESCmkYcJaRNWDeBjJMZvNLzFSLCI7yOIyIiEtQ04igha+c3O5m2ZRoDbxzILdG3eB1HREQk\n6GnEUUJSls1i8IrBVCxdkYntJnodR0TEMzb7a+wzI+BJrWErhaPCUULSvG3z2Jyxmfnd51OpTCWv\n44iIeOaHmOvhczjRpD2Uv8LrOBKoUlNdNVPhKCHnyA9HSFqTRJvYNvRt1NfrOCIinsoq6Wx48OXd\nI6CTHtuRi3BZOOoZRwk5SWuT+D7ze2Z3mo0xxus4IiIiIUOFo4SUDQc3sGD7Aka0GEH9qvW9jiMi\nIhJSdKtaQkbm+Uzi0+KJjYplTFs9AC4iktO+fbCujNcpJNipcJSQMXXzVPYe2cuy3ssoE6HfjiIi\nAKVLO1+nTYNph7zNIsFPhaOEhAPHD5C8Ppnbr7udrnW7eh1HRCRgXH898BHMmgUNorxOI4HqV79y\n106FowQ9ay2PrHyEMBPGjA4zvI4jIhJQfp4j2Lgx3BLtbRYJfiocJegt+WQJaf9O4+n2TxMdpd+K\nIiIi/qJZ1RLUTmWeInF1Ig2rNiSxWaLXcUREREKaRhwlqI1fN56Mkxm80vMVIsIjvI4jIiIS0jTi\nKEFr5zc7mb5lOg82fpBborUbgoiIiL+pcJSglGWzGLxiMBVLV2Riu4lexxERESkWdKtagtLcj+ay\nOWMzC7ov4PLSl3sdR0REpFjQiKMEnSM/HGHU2lG0iW3D/Y3u9zqOiIhIsaHCUYLOyDUj+T7ze1I7\np2J+XqBMRERE/E6FowSV9QfWs3DHQkbeMpJ6Vep5HUdERKRYUeEoQSPzfCbxafHUqFCDJ9o84XUc\nERGRYkeTYyRoTN08lY+PfsyKPisoE1HG6zgiIiLFjkYcJSj857v/kLw+mduvu53OdTp7HUdERKRY\nUuEoAc9aS8KqBMJMGDM6zPA6joiISLGlW9US8JZ8soS0f6fxdPuniY6K9jqOiIhIsaURRwlopzJP\nkbg6kRuq3UBis0Sv44iIiBRrGnGUgDZ+3XgyTmawuOdiIsIjvI4jIiJSrGnEUQLWjq93MH3LdB5s\n/CAtolt4HUdERKTYc1U4GmM6GGP2GWP2G2NG5/J5SWPMK9mfv2eMqeHroFK8ZNks4tPiqVi6IhPb\nTfQ6joiIiOCicDTGhAPPAB2BekAfY8yFW3YMAL6z1tYCpgEpvg4qxcvcj+ayOWMzT7d/mstLX+51\nHBEREcHdM45Ngf3W2s8BjDEvA92BvTnadAfGZ3//GvAXY4yx1tqLnfTgiYMMWj6oQKEl9L2691Xa\nxrbl/kb3ex1FRCQkTHxnIleUu8LrGBLk3BSOVwOHcrzOAJpdrI219idjzAmgEnA0ZyNjzCBgEED4\nVeGs+HRFAWNLqKtevjpzuszBGON1FBGRoFb78trUvrw2W7/c6nUUCQFuCsfc/nJfOJLopg3W2ueA\n5wCaNGlitw7X/4lFRET8qWbFmnya8KnXMSTAmRHuBmrcTI7JAHKuulwd+PJibYwxJYAo4FtXCURE\nREQkKLgpHD8AahtjahpjIoHewLIL2iwD+mZ/3xN4+1LPN4qIiIhI8MnzVnX2M4uPAG8A4cA8a+0e\nY0wysNVauwyYC7xgjNmPM9LY25+hRURERKToudo5xlq7Elh5wXtjc3x/BrjTt9FEREREJJBo5xgR\nERERcUWFo4iIiIi4osJRRERERFxR4SgiIiIirqhwFBERERFXVDiKiIiIiCsqHEVERETEFRWOIiIi\nIuKKCkcRERERcUWFo4iIiIi4osJRRERERFxR4SgiIiIirqhwFBERERFXVDiKiIiIiCsqHEVERETE\nFRWOIiIiIuKKCkcRERERccVYa73p2JgjwEFPOs+fysBRr0OIz+h6hhZdz9Ci6xladD2DS6y1tkpe\njTwrHIOFMWartbaJ1znEN3Q9Q4uuZ2jR9Qwtup6hSbeqRURERMQVFY4iIiIi4ooKx7w953UA8Sld\nz9Ci6xladD1Di65nCNIzjiIiIiLiikYcRURERMQVFY4iIiIi4ooKxwsYYy43xqwxxvw7+2vFS7Qt\nb4z5whjzl6LMKO65uZ7GmDhjzGZjzB5jzE5jzF1eZJWLM8Z0MMbsM8bsN8aMzuXzksaYV7I/f88Y\nU6PoU4pbLq7nMGPM3uyfx7eMMbFe5BR38rqeOdr1NMZYY4yW6AliKhx/aTTwlrW2NvBW9uuLeRJY\nXySppKDcXM/TwP3W2vpAB2C6MaZCEWaUSzDGhAPPAB2BekAfY0y9C5oNAL6z1tYCpgEpRZtS3HJ5\nPbcBTay1NwCvAZOKNqW45fJ6Yoy5DEgE3ivahOJrKhx/qTuwMPv7hUCP3BoZY24CqgFvFlEuKZg8\nr6e19lNr7b+zv/8SOAzkuXq+FJmmwH5r7efW2kzgZZzrmlPO6/wa8BtjjCnCjOJentfTWvsva+3p\n7JdbgOpFnFHcc/PzCc5AyyTgTFGGE99T4fhL1ay1XwFkf616YQNjTBgwBRhZxNkk//K8njkZY5oC\nkcBnRZBN3LkaOJTjdUb2e7m2sdb+BJwAKhVJOskvN9czpwHAKr8mksLI83oaY24Eoq21K4oymPhH\nCa8DeMEYsxa4IpeP/uDyFL8HVlprD2lQw3s+uJ4/n+dK4AWgr7U2yxfZxCdy+yG7cB0xN20kMLi+\nVsaYe4EmQFu/JpLCuOT1zB5omQb0K6pA4l/FsnC01ra72GfGmG+MMVdaa7/KLiQO59KsBdDaGPN7\noBwQaYw5Za291POQ4ic+uJ4YY8oDacAT1totfooqBZMBROd4XR348iJtMowxJYAo4NuiiSf55OZ6\nYoxph/Mff22ttWeLKJvkX17X8zKgAbAue6DlCmCZMaabtXZrkaUUn9Gt6l9aBvTN/r4vsPTCBtba\ne6y1MdbaGsAI4G8qGgNWntfTGBMJ/BPnOr5ahNnEnQ+A2saYmtnXqjfOdc0p53XuCbxttbtBoMrz\nembf2nwW6GatzfU/9iRgXPJ6WmtPWGsrW2trZP/N3IJzXVU0BikVjr80EWhvjPk30D77NcaYJsaY\n5z1NJgXh5nr2AtoA/Ywx27P/xXkTVy6U/cziI8AbwMfAYmvtHmNMsjGmW3azuUAlY8x+YBiXXg1B\nPOTyek7GuZvzavbP44X/oSABwuX1lBCiLQdFRERExBWNOIqIiIiIKyocRURERMQVFY4iIiIi4ooK\nRxERERFxRYWjiIiIiLiiwlFEREREXFHhKCIiIiKuqHAUEcknY8zNxpidxphSxpiyxpg9xpgGXucS\nEfE3LQAuIlIAxpingFJAaSDDWvtnjyOJiPidCkcRkQLI3pf3A+AMcIu19rzHkURE/E63qkVECuZy\nnP2UL8MZeRQRCXkacRQRKQBjzDLgZaAmcKW19hGPI4mI+F0JrwOIiAQbY8z9wE/W2peMMeHAu8aY\nX1tr3/Y6m4iIP2nEUURERERc0TOOIiIiIuKKCkcRERERcUWFo4iIiIi4osJRRERERFxR4SgiIiIi\nrqhwFBERERFXVDiKiIiIiCv/BxOLKkvBFkXfAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fbfdf660160>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "tube = ShockTube(1., 0.1, 1., 0.125)\n",
    "x, p, rho, u, x_min, x_max = tube.get_points(0.2)\n",
    "plt.plot(x, p, linewidth=1.5, color='b', label=r'$p$')\n",
    "plt.plot(x, rho, linewidth=1.5, color='r', label=r'$\\rho$')\n",
    "plt.plot(x, u, linewidth=1.5, color='g', label=r'$u$')\n",
    "plt.xlabel('x')\n",
    "plt.legend()\n",
    "plt.axis([x_min, x_max, -0.1, 1.1])\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}