{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbAAAAEYCAYAAAA9AaOpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHnpJREFUeJzt3XuUXWWd5vHvU6kQAsGEINckcg0O8bIEuhPQ6WUpEgMu\nIS5Bgi3STKZbh/GyZmwFHMdA6ziT1dB0Mg6oA0pkkICiaCtCIlJod8tFATuIMWATmwQImBBuAuby\nmz/et6idopKcSp1d++xdz2ets7LPuy/nrZfiPPXu3z77KCIwMzOrm66qO2BmZrYrHGBmZlZLDjAz\nM6slB5iZmdWSA8zMzGrJAWZmZrXUUQEm6auS1klaUWibLGm5pFWSlkmaVFh3gaQHJa2UNLvQfqyk\nFXndokL7OEnX5fY7JB1cWHd2fo1Vkj5YaD9U0p15n6WSxpY7CmZm1oqOCjDga8CcAW3nA8sj4kjg\n1vwcSTOAM4AZeZ/LJCnvczkwPyKmA9Ml9R1zPrA+t18KLMzHmgx8FpiZHwskTcz7LAQuyfs8lY9h\nZmYV66gAi4ifkkKi6BRgSV5eAszNy6cC10bEpohYDTwEzJJ0ILBXRNyVt/t6YZ/isW4ATsjL7wSW\nRcTGiNgILAdOyoH4NuBbg7y+mZlVqKMCbDv2j4h1eXkdsH9ePghYU9huDTBlkPa1uZ387yMAEbEZ\neFrSPjs41mRgY0RsHeRYZmZWoToE2Msi3fdqpO595XtsmZl1sO6qO9CCdZIOiIjH8+nBJ3L7WmBa\nYbuppJnT2rw8sL1vn9cAj0rqBiZGxHpJa4Gewj7TgB8DG4BJkrryLGxqPsY2JDnszMx2QURo51sN\nrg4zsO8BZ+fls4EbC+3zJO0m6VBgOnBXRDwOPCNpVq5hnQV8d5BjnUa6KARgGTBb0iRJewMnArfk\nGd9twOmDvP42ImLnD9K0bmfbLViwoLXjDWH77W0zWHsrbcXn21tux8Nj4bHwWGz7WLQomDmznLH4\n1KeCE07Y/vqX9y28lw1nLIaro2Zgkq4F3gq8WtIjpCsD/xdwvaT5wGrgfQAR8YCk64EHgM3AudE/\nIucCVwHjgZsi4ubcfiVwtaQHgfXAvHysDZI+B9ydt7so0sUcAOcBSyV9HrgnH6NUPT09bd9+e9sM\n1t5KW/H5UPs7FB6LHfdluNt7LHa+TSeOxV139S+3cyzGjIGDD97++o4bi3b+VTRaH2kYWwDp0WAL\nFiyougsdw2PRz2PRb7hjsWhRxMyZ7enLQJ/5TMRFF7WwYZvey/J75y6/99bhFKLVSJl/ddeNx6Kf\nx6JfO8Zi1qzh92MwY8bA1q07365TOMCsrfxG1c9j0c9j0a8dY3HnncPvx2C6umDLlnKOXQYHmJlZ\nzRRrYO00ZowDzMzMSjRzZjnHdYCZmVmpXANLHGBmZjXjGljiADMzqxnXwBIHmJlZzbgGljjAzMxq\nxjWwxAFmZlYzroElDjAzs5pxDSxxgJmZ1YxrYIkDzMysZsqsgTnAzMysNGXVwHwRh5mZlaqsGpgv\n4jAzs1K5BpY4wMzMasY1sMQBZmZWM66BJQ4wM7OacQ0scYCZmdWMa2CJA8zMrGZcA0scYGZmNeMa\nWOIAMzOrGdfAEgeYmVnNuAaWOMDMzGrGNbDEAWZmVjNl1sAcYGZmVpoya2C+iMPMzErjGljiADMz\nqxnXwBIHmJlZzbgGljjAzMxqpqwamD/IbGZmpSqrBuYPMpuZWalcA0scYGZmNeMaWOIAMzOrGdfA\nEgeYmVnNlPk5sM2byzl2GRxgZmY1U1YNbOxYB5iZmZWorBpYdzds2lTOscvgADMzq5myamCegZmZ\nWanKqoF5BmZmZqUqqwbW3e0ZmJmZlaisGphPIZqZWanK/BzYli0QUc7x280BZmZWM2XVwKR6nUZ0\ngJmZ1UxZNTCo14UcDjAzs5opqwYG9aqDOcDMzGqmrBoYeAZWCkmrJf2LpHsl3ZXbJktaLmmVpGWS\nJhW2v0DSg5JWSppdaD9W0oq8blGhfZyk63L7HZIOLqw7O7/GKkkfHKmf2cxsMGXVwMAzsLIE0BMR\nR0dE33++84HlEXEkcGt+jqQZwBnADGAOcJkk5X0uB+ZHxHRguqQ5uX0+sD63XwoszMeaDHwWmJkf\nC4pBaWY20lwDS+oUYAAa8PwUYEleXgLMzcunAtdGxKaIWA08BMySdCCwV0T0TcC/XtineKwbgBPy\n8juBZRGxMSI2AstJoWhmVokya2C+CrEcAfxI0s8l/WVu2z8i1uXldcD+efkgYE1h3zXAlEHa1+Z2\n8r+PAETEZuBpSfvs4FhmZpUoswZWp1OI3VV3YAjeEhGPSdoXWC5pZXFlRISkmnz8zsxs15VZA6vT\nKcTaBFhEPJb/fVLSd0j1qHWSDoiIx/PpwSfy5muBaYXdp5JmTmvz8sD2vn1eAzwqqRuYGBHrJa0F\negr7TAN+PLB/F1544cvLPT099PT0DNzEzKwtyqyBlTkD6+3tpbe3t23HU9TgniGS9gDGRMSzkvYE\nlgEXAe8gXXixUNL5wKSIOD9fxPENUshNAX4EHJFnaXcCHwPuAn4ALI6ImyWdC7whIv6TpHnA3IiY\nly/i+DlwDKkG9wvgmFwP6+tftDSOfdeR1GDMzawzLV4M11xTXh3s6KPhyivhmGN2sFGb3sskERED\nr21oWV1mYPsD38kXEnYD10TEMkk/B66XNB9YDbwPICIekHQ98ACwGTi3kDDnAlcB44GbIuLm3H4l\ncLWkB4H1wLx8rA2SPgfcnbe7qBheZmYjzTWwpBYBFhEPA28apH0DaRY22D5fAL4wSPsvgDcM0v4S\nOQAHWfc14GtD67WZWTlcA0vqdBWimZlR3xpYuznAzMxqxp8DSxxgZmY143shJg4wM7Oa8b0QEweY\nmVnN+F6IiQPMzKxm/H1giQPMzKxmXANLHGBmZjXjGljiADMzqxnXwBIHmJlZzfhzYIkDzMysZsq+\nF+If/1je8dvJAWZmVjNl1sB2390BZmZmJSmzBjZuHLz4YnnHbycHmJlZzZRZAxs3Dl56qbzjt5MD\nzMysZsqsgTnAzMysNGXXwHwK0czMSlF2DcwzMDMzK4VrYIkDzMysZsqsgfkUopmZlabMGphnYGZm\nVhrXwBIHmJlZzbgGljjAzMxqxjWwxAFmZlYzroElDjAzs5pxDSxxgJmZ1UyZNTCfQjQzs9L4XoiJ\nA8zMrGZcA0scYGZmNVNmDWz33eGFF8o7fjs5wMzMaqbMGtiee8Lzz0NEea/RLg4wM7OaKbMGtttu\nMGZMPU4jOsDMzGqmzBoYwIQJ8Nxz5b5GOzjAzMxqpswaGMBee8Gzz5b7Gu3gADMzq5kya2DQ4BmY\npN0ljSujM2ZmtnNl1sCgPgHWvbMNJHUBc4EzgTeTQk+StgA/A64BboyowzUrZmb1V3YNrEmnEHuB\nY4GLgcMi4sCIOAA4LLf9KXB7aT00M7NtlF0Da8wMDDgxIl5xQWVuuwO4w6cUzcxGjmtgyU4DrBhe\nkiYD04FxhfU/GSzgzMysHGXXwOpyCrGVGRgAkv4S+BgwFbgPOI5UA3t7OV0zM7PBlF0D23tv2LCh\n3Ndoh6FchfhxYCbwu4h4G3A08HQpvTIzs+0quwa2777w5JPlvkY7DCXAXoyIFyBdSh8RK4HXltMt\nMzPbnrJrYHUJsJZPIQJrJO0N3Agsl/QUsLqUXpmZ2XaVXQPbb7+GBVhEzM2LF0rqBV4F3FxGp8zM\nbPvKroHtuy888US5r9EOrXyQeTzwYeAI4F+AKyOit+R+mZnZdpRdA9tvv3oEWCs1sCWkDzKvAE4G\nLim1R2ZmtkNl18AOOACefhr+8IdyX2e4WjmFeFREvAFA0hXA3eV2yczMdqTsGlhXFxx2GDz0ELzx\njeW+1nC0MgPb3LcQEZt3tGFTSZojaaWkByWdV3V/zGx0K7sGBnDkkfCb35T/OsPRSoC9UdKzfQ/g\nDYXnz5TdwapJGgN8EZgDzADOlHRUtb0ys9Gs7BoYwJ/8CfzsZ+W/znDsNMAiYkxE7FV4dBeWXzUS\nnazYTOChiFgdEZuApcCpFffJzEaxsmtgALNnw/e/D1u3lv9au2oonwMbraYAjxSerwFe8ffPqlWD\n7ywVnx2R/nlwsHXb22f460Zqnyb2od3HGzcOxo9PNQazXVV2DQzSDOzVr4bzzoOzzkq/twARsJmj\n2MRYNv8CNm3qf2zenP596SXYuBGeeio9NmyANWvg4YdhwQI47bT29LGVy+j/AQhgsP8lIyJOaU9X\nOlZL33N2/PEXvrw8fnwPe+zRwyu/Ie2m9M/JDLIuv9gOXm1X1o3UPk3sQ7uPF5H+x37xxfRmsOee\n6TFhQrrq65BD4NBDU+3h+ONhypTtv4aNbiNRA5PghhvgU5+C978//e726eYGxrKJsX8F3d0wdmz/\no7sbdtsNJk1K91ScPBlmzEgzut//vpd77unl/vvb1MedfQ+lpCdJs45rgb6Ja1+YRUQ0+rvAJB0H\nXBgRc/LzC4CtEbGwsE1r3+fZ92e5v/tzVNu6NV2e/Pzz6fHcc/DYY+mv04cfhgceSLWHAw6AM8+E\nD30o/SVsBrB4cbo6cPHiCjvRpvcySUTEDs5l7FgrpxAPBE4kfSPzmcAPgGsj4le7+qI183NguqRD\ngEeBM0jjYLZLurrSrGvChP62gZcqb92aQuyrX00zsk9+Ev76r9NfuGYjUQOrg1Yu4tgcET+MiA+S\nvkLlIeB2SR8pvXcdIH904CPALcADwHUR8etqe2VN19UFb3kLXHkl3HMP3HYbvP3t9fiKCyvfSNTA\n6qClUrKk3SW9F/h/wH8GFgHfKbNjnSQH+Gsj4oiI+J9V98dGl0MOgZtvTpdOv/Wt6Q4JNrqNRA2s\nDlq5iONq4HWkKxD+JiJWlN4rM9tGVxdcfHG6AGTevHR585gxVffKqjISnwOrg1ZmYH8OTCd9oeU/\nFz/UPBo+yGzWSf7+79NXvX/5y1X3xKrkGliy0xlYRPgTK2YdorsbvvKVdCrxPe+BAw+sukdWBdfA\nkp2Gk7Sjj2u2vo2ZtceMGfCBD8Df/m3VPbGquAaWtDK76pX0SUlHDlwh6bX55raN/iyYWaf5xCfg\nqqtg/fqqe2JVcA0saSXAZgPrgf8j6TFJq/Jd2R8j3eR2HfCOMjtpZtuaOhXmzoUrrqi6J1YF18CS\nnd6JY5uN053Z++4J8PuI2FJKr2rGd+KwKvzjP8KHPwwrVuz4vozWLIsXw8c/XvHbSIfciWNIF2hE\nxJaIWJcfDi+zCr35zelWVL/8ZdU9sZHmGljiKwzNaqqrK90rcenSqntiI801sMQBZlZj73pXukuH\njS6ugSVDCjBJn5B0q6RfSfqCJN9a1KxCs2bB736X7mZvo4c/B5YMdQb2m4g4AXg9cCvw39vfJTNr\nVXc3nHACLFtWdU9sJLkGlgw1wA6QdDKwZ0TcCtxdQp/MbAje8Y50t3obPVwDS1r5PrCiacAk4BxJ\n+wDdkiYCU4pf8GhmI+f442HRoqp7YSPJNbBkqDOw7wI/i4jTI+LtwDmkb2c+ue09M7OWvO51sGYN\nPPVU1T2xkeIaWDLUz4HdExH/VHj+24i4GpjX9p6ZWUu6u+HYY+Fun9AfNVwDS9pyGX1E+Booswod\ndxzccUfVvbCR4hpY4s+BmTXA0Uf7jhyjiWtgiQPMrAFe/3q4//6qe2EjxTWwxAFm1gDTp8O//Ru8\n8ELVPbGR4BpY4gAza4DddoPDD4eVK6vuiY0E18ASB5hZQ/g04ujhGljiADNriNe9zgE2WrgGljjA\nzBpi+nT47W+r7oWNBNfAEgeYWUMcfrgDbLRwDSxxgJk1xOGHw7/+a8VfNW8jwjWwxAFm1hCTJ6dv\naV6/vuqeWNlcA0scYGYNcthhaRZmzeYaWOIAM2sQ18FGB9fAEgeYWYM4wEYH18ASB5hZg0yblr4b\nzJrNNbDEAWbWIFOnOsBGA9fAEgeYWYNMnQpr11bdCyuba2CJA8ysQaZM8QxsNHANLHGAmTXIvvvC\nM8/Aiy9W3RMrk2tgiQPMrEG6uuCgg3waselcA0scYGYNM2WKA6zpXANLHGBmDeMrEZvPNbDEAWbW\nMAcdBI8+WnUvrEyugSUOMLOG2W8/eOKJqnthZXINLHGAmTXM/vs7wJrONbDEAWbWMPvtB+vWVd0L\nK5NrYIkDzKxhfAqx+VwDSxxgZg3jAGs+18ASB5hZw/QFWETVPbGyuAaWOMDMGmb8eBg3Dp5+uuqe\nWFlcA0scYGYN5NOIzeYaWNLxASbpQklrJN2bHycV1l0g6UFJKyXNLrQfK2lFXreo0D5O0nW5/Q5J\nBxfWnS1pVX58sNB+qKQ78z5LJY0diZ/bbDgcYM3mGljS8QEGBPB3EXF0fvwQQNIM4AxgBjAHuEyS\n8j6XA/MjYjowXdKc3D4fWJ/bLwUW5mNNBj4LzMyPBZIm5n0WApfkfZ7KxzDraA6wZnMNLKlDgAFo\nkLZTgWsjYlNErAYeAmZJOhDYKyL6JtlfB+bm5VOAJXn5BuCEvPxOYFlEbIyIjcBy4KQciG8DvpW3\nW1I4llnH2mcf2LCh6l5YWVwDS+oSYB+V9EtJV0qalNsOAoq3LF0DTBmkfW1uJ//7CEBEbAaelrTP\nDo41GdgYEVsHOZZZx5o82QHWZK6BJd1VdwBA0nLggEFW/TfS6cC/yc8/B1zCyJzGG9JFyBdeeOHL\nyz09PfT09LS5O2atc4A1W11rYL29vfT29rbteB0RYBFxYivbSboC+If8dC0wrbB6KmnmtDYvD2zv\n2+c1wKOSuoGJEbFe0lqgp7DPNODHwAZgkqSuPAubmo/xCsUAM6va5Mnw8MNV98LKUtca2MA/7i+6\n6KJhHa/jTyHmmlaf9wAr8vL3gHmSdpN0KDAduCsiHgeekTQr17DOAr5b2OfsvHwacGteXgbMljRJ\n0t7AicAtERHAbcDpebuzgRvb/kOatZlnYM3mGljSETOwnVgo6U2kU3oPAx8CiIgHJF0PPABsBs7N\ngQNwLnAVMB64KSJuzu1XAldLehBYD8zLx9og6XPA3Xm7i/LFHADnAUslfR64Jx/DrKM5wJrNNbBE\n4fvNDJukaGkc+67y95hbye69F845B+67r+qeWLstXgzXXFPxLKxN72WSiIjBrjJvScefQjSzofMM\nrNnqWgNrNweYWQM5wJrNNbDEAWbWQBMmwB//CC+9VHVPrAyugSUOMLMGktIs7Kmnqu6JlaGunwNr\nNweYWUPtvTesX191L6wMroElDjCzhpo0CZ55pupeWBlcA0scYGYN9apX+Ustm8o1sMQBZtZQEyc6\nwJrKNbDEAWbWUBMn+hRiU7kGljjAzBrKM7Dmcg0scYCZNZRrYM3lGljiADNrKJ9CbC7XwBIHmFlD\n+RRic7kGljjAzBrKpxCbyzWwxAFm1lA+hdhcroElDjCzhvIpxOZyDSxxgJk1lE8hNpdrYIkDzKyh\nPANrLtfAEgeYWUO5BtZcroElDjCzhho3DrZu9ZdaNpFrYIkDzKyhpPTNzM8/X3VPrN1cA0scYGYN\nNmECPPdc1b2wdnMNLHGAmTWYA6yZXANLHGBmDbbXXg6wJnINLHGAmTWYZ2DN5BpY4gAza7AJE+DZ\nZ6vuhbWba2CJA8yswTwDaybXwBIHmFmDOcCayTWwxAFm1mAOsGZyDSxxgJk1mAOsmVwDSxxgZg3m\nAGsm18ASB5hZgznAmsk1sMQBZtZgDrBmcg0scYCZNZgDrJlcA0scYGYN5gBrJtfAEgeYWYM5wJrJ\nNbDEAWbWYL6VVDO5BpY4wMwabI894A9/qLoX1m6ugSUOMLMGc4A1k2tgiQPMrMEcYM3kGljiADNr\nMAdYM7kGljjAzBps7FiQYNOmqnti7eQaWOIAM2s4z8KaxzWwxAFm1nAOsOZxDSxxgJk1nAOseVwD\nSxxgZg3nAGse18CSjggwSadL+pWkLZKOGbDuAkkPSlopaXah/VhJK/K6RYX2cZKuy+13SDq4sO5s\nSavy44OF9kMl3Zn3WSppbGHd4tz+S0lHlzcKZuVwgDWPa2BJRwQYsAJ4D/CTYqOkGcAZwAxgDnCZ\nJOXVlwPzI2I6MF3SnNw+H1if2y8FFuZjTQY+C8zMjwWSJuZ9FgKX5H2eysdA0snAEbn9r/Jr2g70\n9vZW3YWO0Slj0QkB1ilj0QnaMRaugSUdEWARsTIiVg2y6lTg2ojYFBGrgYeAWZIOBPaKiL6/Q74O\nzM3LpwBL8vINwAl5+Z3AsojYGBEbgeXASTkQ3wZ8K2+3pHCsU/uOFRF3ApMk7T/sH7jB/EbVr1PG\nwgHWWdoxFq6BJR0RYDtwELCm8HwNMGWQ9rW5nfzvIwARsRl4WtI+OzjWZGBjRGwd5FgH9R2rsM/U\n4f1IOzfUX/BWtt/eNoO1t9JWfF7mm5PHYsd9aWX7HQXYaBuLXdmmE8eiWANr91jsaH2njcWIBZik\n5blmNfDx7pHqwyCihW004Hkr+wzLaP+fc2d9Ge72o20sHGDN+70o1sBGc4ARER3zAG4Djik8Px84\nv/D8ZmAWcADw60L7mcDlhW2Oy8vdwJN5eR7wpcI+XybV1wQ8CXTl9uOBm/Pyl4B5hX1WAvsP0u/w\nww8//PBj6I/hZEY3nac44/ke8A1Jf0c6rTcduCsiQtIzkmYBdwFnAYsL+5wN3AGcBtya25cBX5A0\nKb/GicB5+Vi3AacD1+V9bywc6yPAUknHkU41rhvY4YgYOEszM7OSKc8gqu2E9B5SAL0aeBq4NyJO\nyus+DfwHYDPw8Yi4JbcfC1wFjAduioiP5fZxwNXA0cB60gxqdV53DvDp/LKfj4gluf1QYCmpHnYP\n8IGI2JTXfZF0BeTzwDkRcU9pA2FmZi3riAAzMzMbqk6/CtHMzGxQDjAzM6slB1iJJP07SZdLul7S\n/Kr7UyVJp0r6Sr5V14lV96dK+dZlV0j6ZtV9qYqkPSUtyb8T76+6P1Xy70O/ob5PuAY2AiR1AUsj\n4n1V96Vq+SrQiyPiP1bdl6pJ+mZEnF51P6og6SxgQ0T8QNLSiJhXdZ+qNpp/HwZq9X3CM7AWSPqq\npHWSVgxon5NvMvygpPO2s++7gR+QrnKsveGMRfYZ4Ivl9nJktGEsGmWI4/HyHXOALSPa0RHg341+\nuzgWrb1PVP3h5To8gD8jXZa/otA2hnRvxkOAscB9wFGkz6RdChw04BjfrfrnqHIsSJ+9WwicUPXP\nUPVYFLb9ZtU/Q4Xj8QHgXXmba6vue5Vj0dTfh138vRjS+4RnYC2IiJ+S7lJfNBN4KCJWR/rM2FLg\n1Ii4OiL+S0Q8KumtkhZJ+jLpLiO1t6tjAXyUdGPl0yR9aGR7XY5h/F5MlvQl4E1N+it8KOMBfBt4\nr6TLSDcMaJShjEVTfx/6DPH34iMM4X2iE+/EURfFUyCQbvS7zT2iI+J24PaR7FRFWhmLxfTfLaXJ\nWhmLDcCHR7JTFRp0PCLiD6QbFIwm2xuL0fT70Gd7Y/FR4H+3ehDPwHadr37p57Ho57HYlsejn8ei\nX1vGwgG269YC0wrPp7Ht17WMJh6Lfh6LbXk8+nks+rVlLBxgu+7npG+CPkTSbqQ72zfuXH6LPBb9\nPBbb8nj081j0a8tYOMBaIOla4J+BIyU9IumcSF+W+RHgFuAB4LqI+HWV/RwJHot+HotteTz6eSz6\nlTkW/iCzmZnVkmdgZmZWSw4wMzOrJQeYmZnVkgPMzMxqyQFmZma15AAzM7NacoCZmVktOcDMzKyW\nHGBmZlZLDjCzGpE0TtLtktSm4z1XWP6SpOMHthde9yeS/J5hHcO/jGb18ufA96N994ArHmcWcMcg\n7UTES8BPgbltel2zYXOAmdXLmcB3AfKdvFdK+pqk30i6RtJsSf8kaZWkP83b/VdJK/Lj44MdVNJR\nwKqdBOP38uubdQR/I7NZh8lfpf4kcFhEXFxoHwO8PiJWFTY/HHgv6Y7edwNnRMRbJJ0CfFrS54C/\nIH2Fexdwp6TbI+K+AS97EvDDnXTtPuDNu/6TmbWXZ2BmHUTS2cCjEfFtYPaA1a8Gnh3Q9nBE/CrP\nnH4F/Ci33w8cAvx74NsR8UJEPA98G/izQV56NnDzjvqWTyN2Sdp9CD+SWWk8AzPrLGcCJ+eLJfYZ\nZP3AizdeKixvBf5YWO77/7u4jxhQ35I0HpgUEY+30L9X7G9WFc/AzDqEpD1J39G3FXg3udZV8Htg\nwhAP+1NgrqTx+fhzc1vR24DbWujfOGBLnomZVc4zMLPOMQt4QdKpwFHA/yiujIgtku6X9NqI+E1f\n84BjxLa7xL2SrgLuym3/NyJ+OWCfk4BvDmjbQ9IjheeXkK5Q/NmQfiKzEvkbmc06hKTPALdHxMAZ\nUnGbvwD2j4iFbXzdXwAzI2LLTrb7AnB3RHynXa9tNhw+hWjWOQ6l/3NY2/MN4F3t+iAzQEQc20J4\njSNdEHJju17XbLg8AzMzs1ryDMzMzGrJAWZmZrXkADMzs1pygJmZWS05wMzMrJYcYGZmVksOMDMz\nqyUHmJmZ1dL/B+t1HMuN709zAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x108e93d50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import CoolProp, numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "AS = CoolProp.AbstractState('HEOS','Water')\n",
    "\n",
    "T = 300\n",
    "rhoc = AS.rhomolar_critical()\n",
    "\n",
    "# Saturated liquid density\n",
    "AS.update(CoolProp.QT_INPUTS, 0, T)\n",
    "rhoL = AS.rhomolar()\n",
    "plt.axvline(rhoL/1000.0,lw=2,c='r')\n",
    "# Saturated vapor density\n",
    "AS.update(CoolProp.QT_INPUTS, 1, T)\n",
    "rhoV = AS.rhomolar()\n",
    "plt.axvline(rhoV/1000.0,lw=2,c='r')\n",
    "\n",
    "AS.specify_phase(CoolProp.iphase_liquid) # Something homogeneous to avoid flash call\n",
    "\n",
    "plt.axvline(rhoc/1000.0,dashes = [2,2])\n",
    "p = []\n",
    "rho = np.linspace(1e0, rhoc*5, 10000)\n",
    "for D in rho:\n",
    "    AS.update(CoolProp.DmolarT_INPUTS, D, T)\n",
    "    p.append(AS.p())\n",
    "    \n",
    "plt.plot(rho/1000.0,np.array(p)/1e6)\n",
    "#plt.yscale('log')\n",
    "plt.ylim(-1e6,10e5) # Undo this to see how crazy it is\n",
    "plt.xscale('log')\n",
    "plt.xlabel(r'$\\rho$ (mol/L)')\n",
    "lab = plt.ylabel(r'$p$ (MPa)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
