{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Miniconda\\lib\\site-packages\\ipykernel\\__main__.py:20: RuntimeWarning: invalid value encountered in double_scalars\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAESCAYAAADwnNLKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4XNWZ5/HvK1m75UVe5d3Y2IDZHMwSICAWA4FA6Axp\nAiQs7gQmdEKWCel0Og1Kz0zI1jPpEALpELawTkMAA3ZjIFYwSwwG78ZrW8K7jGzJtiRrPfPHUUll\nWbarVKU6Jev3eZ7z3FtVt+59fW3ft+7ZrjnnEBGRvisjdAAiIhKWEoGISB+nRCAi0scpEYiI9HFK\nBCIifZwSgYhIHxckEZjZd8xshZktN7MnzSwnRBwiIhIgEZjZaOCbwGnOuZOATOBLqY5DRES8fgGP\nm29mLUA+sCVQHCIifV7K7wicc1uAfwU+BrYC1c6511Mdh4iIeCGqhgYDVwETgFFAfzO7IdVxiIiI\nF6Jq6GJgo3OuCsDM/gScDTwR2cDMNAGSiEg3OOcs3u+E6DVUAZxlZnlmZvjEsKrzRs65Xlvuvvvu\n4DEo/vBx9LXYFX9Huewyx8svpz7+7grRRvAe8CzwIbCs7e1/T3UcIiI9pbwcJkwIHUXsgvQacs6V\nAqUhji0i0pOcg4oKGD8+dCSx08jiHlBSUhI6hIQo/nB6c+yg+AF27oT8fOjfP/F4UsUSqVfqKWbm\n0jEuEZEjee89uP12WLQo9cc2M1wvaSwWETlq9bb2AVAiEBFJKiUCEZE+rrc1FIMSgYhIUumOQESk\njwuVCO5ZcE+3v6tEICKSJM75RBCiaujVDa92+7tKBCIiSVJVBTk5MGBA6o9dXl3e7e8qEYiIJEmo\naqHm1ma27dvW7e8rEYiIJEmoRLCpZhMjCkZ0+/tKBCIiSRIqEZRXlzNx8MRuf1+JQEQkSUKNISiv\nLmfCoAnd/r4SgYhIkoS6I9hYvZEJA7t/YCUCEZEkUdWQiEgfFnIMgaqGRETSwO7dkJkJgwal/tgb\nqzcqEYiIhBaqWqixpZHK2krGDBjT7X0oEYiIJEHIMQSjCkfRL6P7Tx5OeSIws6lmtjiq1JjZHamO\nQ0QkmUI2FCdSLQQBHl7vnFsDTAcwswxgC/B8quMQEUmmioqAXUcTTAShq4YuBjY45zYFjkNEJCFB\nu44O6n7XUQifCL4EPBk4BhGRhPXmqqFgicDMsoErgf8IFYOISLIEHVXc29oIonwW+MA5t7OrD0tL\nS9vXS0pKKCkpSU1UIiJxqq72A8pSPYagrKyMZc8s4/m1z/Pn3D93ez/mnEtiWHEc2OxpYK5z7tEu\nPnOh4hIRideSJXDjjbBsWWqP29DcwICfDqDuh3VkZmRiZjjnLN79BKkaMrMCfEPxn0IcX0QkmUJV\nC1XUVDBmwBgyMzIT2k+QqiHnXC0wNMSxRUSSrTf3GILwvYZERHq93vocggglAhGRBPXmrqOgRCAi\nkrCQXUdVNSQikgZ0RyAi0ofV1EBTExQVpf7YSgQiImkgMtmcxd17PzH1TfXsrt9NcWFxwvtSIhAR\nSUDIMQTjBo4jwxK/jCsRiIgkoLe3D4ASgYhIQkKNIdi4O/HJ5iKUCEREEtDbRxWDEoGISEI2bgyU\nCGpUNSQiElxrK6xdC1Onpv7YqhoSEUkD5eV+/MCAAQGOXV3OxMGqGhIRCWrlSpg2LfXHrW2sZW/j\nXkYUjEjK/pQIRES6KVQiqKipYPzA8ViSRrEpEYiIdFOoRJDM9gFQIhAR6bZQiSCZXUdBiUBEpFta\nWmDNGjjhhNQfO5mjikGJQESkWzZuhGHDoH//AMeuPgqqhsxskJk9a2YfmdkqMzsrRBwiIt0VqloI\nktt1FAI9vB74N2COc+4aM+sHFASKQ0SkW0Ingl59R2BmA4HPOOceAnDONTvnalIdh4hIIkIlgr0N\ne6lvrmdY/rCk7TNE1dBEYKeZPWxmH5rZ780sP0AcIiLdFrLH0IRBE5I2hgDCVA31Az4FfMM5976Z\n/Qr4AXBX9EalpaXt6yUlJZSUlKQwRBGRQ2tp8XMMHXdc6o8dXS1UVlZGWVlZwvs051zCO4nrgGYj\ngXedcxPbXp8L/MA597mobVyq4xIRidXatXDppb7nUKrdu/BeVn+ymvuuuO+gz8wM51zctwoprxpy\nzm0HNpnZlLa3LgZWpjoOEZHuCtlQvLF6Y1J7DEG4XkPfBJ4ws2xgA3BLoDhEROIWusfQ2WPPTuo+\ngyQC59xS4PQQxxYRSdTKlfDZz4Y5drK7joJGFouIxC101ZASgYhIQM3NsG4dHH986o9dvb+a5tZm\nhuQNSep+lQhEROKwfj2MGgX5AUY/VVRXJH0MASgRiIjE5WirFgIlAhGRuITuMZTM5xBEKBGIiMQh\nZCL4aOdHTBky5cgbxkmJQEQkDiETwfLK5Zw84uSk71eJQEQkRk1NsGFDmDmGWl0ryyuXc9Lwk5K+\nbyUCEZEYrVsHY8ZAXl7qj11eXc6g3EEMzhuc9H0rEYiIxChotdCOnqkWAiUCEZGYhUwEy3Ys4+Th\nSgQiIkEFTQSVyzhpRPLbB0CJQEQkZqoaEhHpwxob4b/+C6ZOTf2x65rqqKipYOqQnjm4EoGISAzW\nroXx4yE3N/XHXrVzFVOGTCErM6tH9q9EICISg6O1WgiUCEREYhK6x1BPDCSLUCIQEYlB6B5DuiMQ\nEQksVCJwzvkxBD2YCII8s9jMyoE9QAvQ5Jw7I0QcIiKxaGiA8nKYkvyJP49oR+0OWl0rxf2Le+wY\nQRIB4IAS59yuQMcXEYnZmjUwcSLk5KT+2JG7gWQ/lSxayKqhnvtTiYgkUfAeQz00tUREqETggNfN\nbJGZfS1QDCIiMQndUNxTU0tEhKoaOsc5t83MhgGvmdlq59yC6A1KS0vb10tKSigpKUlthCIibZYu\nhZtuCnPs5TuW8/en/32Xn5WVlVFWVpbwMcw5l/BOEgrA7G5gn3PuX6Pec6HjEhEBcA6GDYMlS/yz\nCFKpubWZAfcMYOedOynILjji9maGcy7uaveUVw2ZWb6ZFbatFwCXAMtTHYeISCzWrYP8/NQnAYC1\nVWsZPWB0TEkgESGqhkYAz7e1gPcDnnDOzQsQh4jIEb3zDpx9dphjxzW1RH19t4+T8kTgnNsInJrq\n44qIdEfIRBDXw2huvLHbx9HIYhGRwwiaCOLpMbR0abePo0QgInII1dV+RPEpp4Q5fsxVQ/v2webN\n3T6OEoGIyCH89a8wYwZk9cxjAA6rZn8Nn9R9wjGDjznyxsuXJzTQQYlAROQQgjYUVy5n2vBpZFgM\nl+klSxK6bVEiEBE5hHffDdxjKNaG4qVL4dTu98FRIhAR6UJLCyxcCGedFeb4y3bE0VCsOwIRkeRb\nsQJGjYKhQ8McP+aH0bS0+GBP7v7EdEoEIiJdCNk+4JxjReWK2B5PuX49DB8OAwd2+3hKBCIiXQiZ\nCCpqKuif3Z8h+UOOvHGC7QOgRCAi0qXgI4pjnVoiwfYBUCIQETnI9u2waxccd1yY46eyxxB0MxGY\nWYGZZSZ0ZBGRNPXuu/DpT0NGoJ/KcU0tkao7AjPLNLPrzewVM6sE1gDbzewjM/uFmU1OKAoRkTQS\nsloI4qga2rkTamth/PiEjhdrvvszMBn4R6DYOTfGOTcMOBdYCPzMzL6SUCQiImkiZCLY37yf8upy\njhsaQ73U0qX+biDBB9vHOg31TOdcY+c3nXNVwLPAs2YWYDYOEZHkamjwtS1nnBHm+Kt2rmJy0WSy\nM7OPvHEkESQopjuCrpKAmRV32qYp4WhERAL78EOYOhX69w9z/EVbFzF95PTYNl6yJOGGYkis19Ar\nCR9dRCTNhG4fKCsvo2RCSWwbp/KO4BASq5QSEUlDoUcUx5wIGhr8A5UTmH46ItZeQ11d9B/stI3G\nJIhIr+Zc2ESwbtc6+mX0Y+KgiUfeeNUqmDQJcnMTPm6sF+/5ZvZNMxsXecM5d5+ZZZvZRWb2GHBT\nPAdu65K62Mxeiud7IiI9pbzcd8BJsDdmt0XuBrr+7d1JktoHIPZeQ58FZgFPmdkxQDWQC2QC84D/\n65xbHOexvwWsAgrj/J6ISI+I3A0k2Buz28rKy7j4mItj2zhJ7QMQe6+heufcfc65c4DxwEXAp5xz\n45xzX403CZjZGOByfPWS2hpEJC30mvYBSOodQdz1+s65RufcVufc7gSO+3+BO4HWBPYhIpJUvaZ9\nwLmk3hHEWjWUNGb2OaDSObfYzEoOtV1paWn7eklJCSUlh9xURCRhe/fC2rUwPcYu/MkWV/vAxx9D\nXh5lq1ZR9tvfJnxsc84lvJO4Dmj2E+ArQDO+nWEA8Jxz7saobVyq4xKRvu2NN+Duu+Gtt8Ic/7rn\nrmPmMTOZNX3WkTeePRvuvx/mzj3gbTPDORd3dXvKu3w6537onBvrnJsIfAn4c3QSEBEJYcGCvtk+\nAOnxPAL99BeR4F5+GS6/PMyx11atJSsjK7b2AUhq+wAkmAjMrH/bMqs7zydwzv3FOXdVIjGIiCRq\n82bYuBHOOSfM8cvKy7hg4gWxtQ9AUh5GE63bicDMvg/cZWb/CgwEHkhaVCIiKfTyy/DZz0JWoDmU\nyyrKKBlfEtvGe/bAtm1w7LFJO34idwQLgbuA7+PHFaRDNZOISNxmz4arAtVNxN0+sHy5n18oM3kP\niTzixdvMMszs62Z2rpllt703GKgFbnbOtTjnngHeSFpUIiIpsm+fbyi+9NIwx19btZbszGwmDJoQ\n2xeS3FAMMSQC51wrUANcDVzR9vb/AkYTNRW1c+7JpEYmIpIC8+b55xMPHBjm+HGNH4CkNxRD7NU5\nmc657znnnm973Qp8BphtZn+b1IhERFIoZLUQxNk+AGHuCNoM6vT6aefc94DTAD2iUkR6pZYWeOUV\nuPLKMMePu32guRlWroSTY3iwfRxiTQRDzaw9GTjn3m5btgIFSY1IRCRF3n0XRo8ON+103O0D69ZB\ncTEUJnfS5lgTwf3Ak2Z2XvSbbQ+jOTGpEYmIpEjwaqE0aB+AGCedc85tN7PbgcfNbABQBjQAZwP/\nlvSoRERSYPZsePzxcMefXz6fSyfF0V1p4UI47bSkxxFz33/nXLlz7lzgNqAC2ALMcs79v6RHJSLS\nw9as8TOOfupTYY4fd/sAwGuvwcUxPrgmDnFPQ+2cexd4N+mRiIik0Esv+UbijEBDYddUrSG3Xy4T\nB8c4v9DWrb6EvCMQETmapEv7QMxeew0uuiipI4ojlAhEpM/55BPf7nrhheFi6FYimDmzR2JRIhCR\nPmfOHP/jOjc3zPHjbh9obVUiEBFJptDVQisqV5CXlRf7+IHly2HAAJgYY3tCnJQIRKRP2b/f/7i+\n4oojb9tTnln5DNccf03sX5g3Dy65pMfiUSIQkT6lrAxOOgmGDQtzfOccT614iutOui72L/VgtRAo\nEYhIHxO6Wui9Le+RlZHF9JHTY/tCfb2fC+OCC3osprjHEYiI9FbO+fEDr70WLoanVjzFdSdeF/u0\nEm+95SeZ68F5slN+R2BmuWa20MyWmNkKMytNdQwi0jctXgx5eTB1apjjt7S28MzKZ+KvFurB9gEI\nkAicc/uBC5xzpwKnApeZ2ZmpjkNE+p5HHoFrr4VYf4wnW1l5GaMLRzNlyJTYvzRvXo+2D0CgqiHn\nXF3bajb+eQatIeIQkb6jthaeeMLfFYQSqRaK2Y4dUF4OZ5zRYzFBoMbitucgLwF2APOcc++HiENE\n+o6nn4ZzzoFx48Icv6G5gedXP8+1J14b+5def903Evfr2d/soe4IWoFTzWwg8LyZTXPOrYzeprS0\ntH29pKSEkpKSlMYoIkeXBx6AH/843PH/c/1/ctLwkxgzYEzsXzpCt9GysjLKysoSjs2ccwnvJKEA\nzP4ZqHPO/WvUey50XCJy9Fi0CL74RVi/vkfmbIvJtc9ey4UTLuS2GbfF9gXn/OPT3nwTJk+O6Stm\nhnMu7haQEL2G2h97aWZ5wEzgo1THISJ9xwMPwK23hksC+xr38er6V7nmhDhGE69aBTk5MGlSzwXW\nJkTVUDHwqJll4hPRM865OQHiEJE+oLoannsOVq8OF8OLq1/k3HHnMiR/SOxfinQbTUEXp5QnAufc\nciDQM4FEpK/54x/h0kthxIhwMTy54kmuP/H6+L40bx7MmtUzAXUSvI2gK2ojEJFkcA6mTYP774fz\nzw8TQ1VdFZN+PYnN391M/+z+sX2pocFPhlReDkVFMR+r17QRiIikyoIFfnneeeFieHbVs1w2+bLY\nkwDAO+/A8cfHlQQSoUQgIket+++H//7fw40kBl8tFNcgMkjJtBLRlAhE5KhUWQn/+Z9w443hYti8\nZzMrKldw2eTL4vtiCqaViKZEICJHpYcfhi98AQYNChfDMyue4W+O+xty+uXE/qVPPoF16+Css3ou\nsE6UCETkqNPaCr/7na8WCqlb1UJz5viW7ezsngmqC0oEInLUmTfPt7Oefnq4GN7Z9A676nfF/oD6\niIceSnl9lhKBiBx1Io3EIf3vBf+bH5zzAzIz4hjOvHYtfPRRyh+hpkQgIkeVjRv9Q72ui7NGJpk+\n3PYhS7cv5eZTb47viw8+CDfdlNJqIdCjKkXkKPOjH8Edd0BBQbgYfrLgJ3zv7O/F10jc2AiPPtox\n+CGFlAhE5Kjx4Ycwf75vKA5l1c5VLPh4AY9e/Wh8X3zpJTjuOJgSx9PLkkRVQyJyVHAOvv99uOsu\n6B/HIN5ku+ete/j2md+mIDvOW5Lf/x6+9rWeCeoINNeQiBwVXn0VvvUtWL4csrLCxLBh1wbOfPBM\nNtyxgYG5A2P/Ynk5zJgBmzZBXl63j6+5hkSkz2pp8XcD99wTLgkA/Oztn3H76bfHlwTAdxm9/vqE\nkkAi1EYgIr3e44/76qCrrw4Xw6aaTTy76lnWfXNdfF9sbvaJYE64x7IoEYhIr1ZfD//8z/7h9CEn\nl/vlO7/k76b/XXwPnwFfpzV6NJx8cs8EFgMlAhHp1e6911evn312uBh27NvBH5f9kZW3r4z/ywEb\niSPUWCwivVZVle9x+dZbMHVquDj+4bV/oLaplt9c/pv4vrhtG5xwgm8kTkJXp+42FuuOQER6rZ/8\nBK65JmwS2FW/iwcXP8ji2xbH/+WHH4YvfjFsf1cCJAIzGws8BgwHHPDvzrlfpzoOEendysvhkUdg\nZTdqY5Lp1wt/zdVTr2bcwHHxfbG1Ff7wB3jqqZ4JLA4h7giagO8455aYWX/gAzN7zTn3UYBYRKSX\n+tGP4JvfhJEjw8VQUV3Bfe/fx7t/9278X54/398JhJwitU3KE4FzbjuwvW19n5l9BIwClAhEJCZz\n58Kbb/pZRkNxzvG1l77Gd8/6LpOLJse/g0gjcciuTm2CDigzswnAdGBhyDhEpPfYvh1mzYI//hEK\nC8PF8dDih9hVv4s7z7kz/i/v2OGfo3nDDckPrBuCNRa3VQs9C3zLObev8+elpaXt6yUlJZSUlKQs\nNhFJT62t/pktX/uaf4hXKJv3bOYHb/yAP9/4Z/pldOMyWloKt9wCgwcnFEdZWRllZWUJ7QMCdR81\nsyzgZWCuc+5XXXyu7qMicpBf/AJeeAH+8hfoF+hnrHOOK568grPGnMVd598V/w5WrIALL4Q1axJO\nBJ31mu6jZmbAH4BVXSUBEZGuvP++TwTvvx8uCQA8tvQxtu7dyj+e+4/xf9k5+B//w7d0JzkJJCJE\nG8E5wJeBC8xscVu5LEAcItJL7N3rnzh2330wfny4OLbu3cqdr93Jw59/mKzMbsxuN3cuVFTA17+e\n/OASoJHFIpL2brwRcnJ8R5tQnHN8/unPc+rIU/mXC/4l/h00Nfn5hH7xC/jc55IfIL2oakhEJB6P\nP+6rgxYtChvHk8ufZGP1Rp7922e7t4Pf/Q7GjIErrkhuYEmgOwIRSVsbNsBZZ8Hrr8Mpp4SLY/u+\n7ZzywCm8cv0rzBg1I/4d7N7tJ0V6/XU46aTkB9hGD6YRkaNKTY2fhudHPwqbBJpamrjlxVuYdeqs\n7iUBgP/5P/3DEnowCSRCdwQiknZqa+HSS32V+n33hRt865zjq7O/yrZ923jxSy92r4F43Tr49Kf9\npEgjRiQ/yChqIxCRo0J9PVx1FUyZAr/5TdgZGO6afxfLKpcx/6b53UsCAHfeCd/7Xo8ngUQoEYhI\n2mhogP/232D4cN9DKCNg5fUDix7g6ZVP8/ast+mf3c1poufPh6VL/ePT0pgSgYikheZmP1YgJwce\newwyM8PF8sLqF/iXv/wLb816i+EFw7u3k8ZG+M534Gc/g9zc5AaYZEoEIhJcSwvcdJOvFnrhBcjq\nZi1MMrz98dvc+tKtzL1hLscMPqb7O/r2t2HcON/ineaUCEQkqNZWuO02/9TGV17xdwShrNq5ii/8\nvy/w+Bce57RRp3V/R7//va8WWrgwLaaZPhIlAhEJprERbr8dVq/2szLn5YWLZcueLVz+xOX8cuYv\nuWTSJd3f0TvvwD/9EyxYAAMGJC/AHqRxBCISxJYtfirpTz6BOXPCPrZ36falnPPQOdxx5h185ZSv\ndH9HW7f6qqCHHw77IOU4KRGISMq9+aZ/QuOVV8Kf/hT2h/OLq1/k4j9ezM8u/hnf/fR3u7+jhgb4\nwhf8LU4aTiNxOBpQJiIp4xz827/BPff4nkGXXhoyFsfP3/459753L89f+zynj07g2cHOwVe/6odD\n/8d/BGsX0IAyEUlrtbX+yWKrV8Nf/woTJ4aLpaG5gVtfvpUVlStY+NWFjB4wOrEd3n+/bxj+6197\nReNwZ6oaEpEet2aNn2UhOxvefjtsEqisreTCxy6krqmON29+M/EksGAB/PjHvt9ryIaOBCgRiEiP\nqa2FH/4QzjnHV50//HDYnkELNy/kjN+fwUUTL+KZa56hILsgsR2+9RZcc42v55o8OTlBBqCqIRFJ\nOud8I/B3v+uTwLJlMGpUuHj2NOzhn974J5796Fnu/ey9XHPCNYnv9E9/8gMgHn88bGNHEigRiEhS\nrVkD3/ym70n56KNQUhI2nhdXv8g35n6DS465hJW3r6Qoryjxnd57L/z0p/Dqq/CpTyW+v8CCJAIz\newi4Aqh0zqXnBN0iEpc9e3xvoN//3o+n+sY3wk4VsXXvVu6YewfLdizjsasf44KJFyS+09ZW+MEP\nYPZs39gxYULi+0wDodoIHgb0wHqRo8CWLfD97/sG4K1bYflyP9daqCTQ0trCA4se4JQHTuH4ocez\n7OvLkpMEGhrgy1/2CeAoSgIQ6I7AObfAzCaEOLaIJMeKFfDLX/ofx1/5CnzwQdhrY1NLE48ve5x7\n3rqH4sJi5t80nxOHn5icndfUwN/8DQwa5B83GbLFuweojUBEYuYclJXBL34Bixf76p/166EoCdXu\n3dXQ3MAjSx7hp2//lGMGH8O/X/nvnD/+fCxZ/fnnz4dbb4XLLoNf/Srs/Ng9RIlARI5oxQp46ilf\ncnN91c+f/hR2mv26pjoe/PBBfv72zzl5xMk88YUnOHvs2ck7wO7dvs7r1Vf98zKvvDJ5+04zaTuO\nwKyUkpJSSktLufnmMszg4ov904vOOssP3vvbv4X33/dPgTOD0lL/3dJSvdZrvU709be/7Rt/R4zw\nz1yfPx+ee87PqXbrrb7TTIj4br+zEjMYdsV9zC+fz+Vb3mPul+cw7w9nJ/d44x6CnBxKb1iHXXVl\n8L+Prl6XlZVRUlLafr3srmBzDbW1EbzUVa+hQ8019NFH/pfJxx8fXPbtg9GjYcwYGDvWL6PL6NH+\n8XchH30nks6am+G993wV+Ny5vsrnmmv8U8POPTfs/53q/dU8veJpHl7yMJv3bObGk2/klum3MGXI\nlOQeaMsWX9+1ejU8+KAfBNGLdHeuoSCJwMyeAs4HhgCVwF3OuYejPu8yEcxZN4eczBwumHgBGXbg\nv8q6Ov93uGkTbN7cUSKvt271d3ojRviBLaNH++WoUVBcDCNH+mVxMQwbpoQhRz/nYN06eO01X8rK\nfGPvzJm+XHBB2O6fjS2NlJWX8ejSR3ll7SvMnDSTWafOYuakmfTL6Jfcg+3bB7/7nb/Nuf12Pxw6\n5BNyuqlXJYIjOVQieGr5U/z8nZ+zq34XN558IzedehOTiybHvN/GRti+3SeMrVs7ltu2+bJ9u19W\nV/tkUFzsE0fnMnKkv7sYPhyGDDkq247kKFRdDYsW+erU99/3v/7N/EX/4ovhoov8v++QdtXvYs66\nOby09iXmbZjH1CFTueGkG7j+pOsZkj8k+QesqvKDw+67z2e+u++GadOSf5wU6ROJgDfegOJilhQ1\n8uiSR3lyxZMcW3QsN596M1884YsMzB2YlOM3NkJlpU8KO3Ycuuzc6e8yBg/2iWP4cL8cNgyGDj24\nDBniS0FBr5ygUHqJlhaoqPBVqR99BB9+6C/827fD9On+OQCRcswxYf8tOudY/clq5qybw+y1s1m8\nbTEXTryQK6dcyRVTrmBk/5E9c+DNm+H//B945BHf8Pj978Oxx/bMsVKobySCRx7xf2FXXw0//jFN\nw4cyd/1cHlnyCG9sfIMLJlzA56d+niumXMHwguEpibWlxf+oqKz0iaGy0j9xKVKqqjrWd+6EXbt8\nXWxRkU8KRUUd64MHH1yKivxy0CAYONDP3ijS0uJ/qFRU+LJuXceFf906/8Pj+ON9OeUUf9E//vjw\nd6+trpXlO5bzZsWb/KXiL7xZ8Sb5WflcOulSrpp6FRdOvJC8rB7qo++c7/P629/6Lk+33OInQxqd\n4OyjaaRvJALwP8F/8hN46CE/ocn3vgf9+1NVV9X+q+K1Da8xbfg0rppyFVdNvYrjhh6XvD7FSbB/\nv08IVVUdy6oq/0frXHbt8suaGl+ysnxSiJSBA/3TnSLLzuuFhV2XkHW/cni1tf7Xe+eyZQuUl/sL\n/+bN/kfC+PG+TJrUceE/7rj0mQ15V/0uPtz2IYu2LuKdTe/w1sdvMTR/KOePP5/zJ5zPeePPY9zA\ncT0bxH/9Fzz5JDzxhB8dfPPN8Pd/7399HWX6TiKIKC/3E5rMn+/7Uc2aBf18A1JDcwNl5WXMXjOb\n2Wtnk9sPuJEkAAAO5ElEQVQvl4smXkTJhBLOH38+xYXFPf5n6AnO+Ubx6mqfFKqrfdmzp6PU1By4\nvndv1yUz018sIqWgoGM9P9+/jiyj1/PyDl1ycw8sOTl9r9G9tdUn+ro6f5737Dl4uWePT/CdS+QH\nQXNzRweG6FJc7C/6EybAuHFh+/B3ZWftTpbtWMairYv4YNsHfLDtAyprK5k+cjozRs3gzNFnct74\n81Lz/2/nTnjmGX/x37DB9zW/4YaOvudHqb6XCCIWLYI77/SV9nfeCdde669abZxzLN2xlLLyMsrK\ny3iz4k2GFwynZEIJJRNKOG/8eYwqDDg/bgDO+R9G+/b5Ulvbsb5vn7+I1dYevKythfr6Q5f9+31p\naOhYz872CSEnx69HXkeWWVl+PSvr4NKvny+ZmQeuR0pGxoHFrGMZ0fn/fGur//O3tnYU53xVS3Nz\n16Wx0f+ZOi8jf866Ol8i5yA31yfGwsKOu7Lo5YABHVWC0VWDkWVhYXpfq6rqqli5cyUrK1f6Zdt6\nY0sjJ404iRnFMzht1GnMGDWDY4uOJTMjBfVRzc2+IWTePN8Favly+Nzn/MV/5sw+cwvcdxMB+P/J\nr74Kv/kNvPuu/8u/7bYuW/9bWltYXrm8PTEs+HgBef3yOH306cwonsHpo0/ntOLTeqaHQh/jnL9o\n7t9/4EW08wW1qanrErkQRy7S0cvoC3nnEn38zvEcKnlkZByYfKJLVtaBySx6PS/P/+7Iz++4K+rt\nd0HOOarqq9iwawPrdq1j/a717cv1u9bT3NrMtGHTfBnesSzuX5y6KljnfJXP66/7i/+f/+xvky65\nxF/4P/OZo24+oFj07UQQraLCDwT5wx98xeltt/lRMYe4j3bOsbF6I+9veZ9FWxexaNsiPtj6AUPz\nh3LaqNMO+Ad/bNGxZGX2jV8WcvRqamli275tbNmzhYqaCiqqK6ioqaC8urz9db+MfkwumsyxQ45l\n8uC2ZdFkJhdNZlj+sNS3uVVX+1/8Cxf68t57/rbwoov8xf/ii33dWR+nRNBZUxO8/DI88ICfFvHy\ny/1cIZde6u/ND6PVtbK2ai2Lty0+4NZ3055NHDP4GKYNm8ZxQ49jctFkJg2exKSiSYwoGJFWDdLS\n9zQ0N1BZW8n2fdvby47aHWzZs4Ute9vKni1U1VcxvGA4owtHM37QeMYP9GXCoAntr5PVFTtuzvmW\n8BUrfFm+3F/0N2/2D4A588yOMmZMetehBaBEcDgffwwvveTL22/7BqMrr/Qljqdo1zfVs6ZqDSsr\nV7Kmag0bdm9gw64NbNi9gfqmeiYVTWLS4EmMHziesQPHMnbAWMYOHMuYAWMo7l+cmrpSOSq0ulb2\nNOyhqq6KXfW7qKqvoqquip11O9lZu5NP6j7x622vK2sr2de4j+EFwxnZfyQj+o9gZIFfjiocxejC\n0YweMJrRhaMZ0X9E8kfmxquhwXf4WL++o+9r5OKfnw8nnujLtGkwY4Zf7xc45l5AiSBWe/f6xqSX\nXoJXXvEdrs89F84+Gz79aZgypVu/Mmr217QnhoqaCjbv2cymPZvYVLOJTXs2UVVXxcj+IykuLGZE\nwQj/n7VghP8P27Y+JH8IQ/KGUJRXpCqoXq6ltYU9DXvaS01DTcf6/hqq91e3l937dx+wvqt+F9X7\nqynIKqAor4iivCKG5Pt/F8Pyh/lS4JdD84cyrGAYwwuGU5RXdNDUK8G0tPg+r5HJwCoqYONGf+Ff\nv94P6R83zj/wPdL3NXLhHzo0dPS9lhJBd7S2+gEm77zTUWprfUI4+2w44ww/7eLwxAenNbY0snXv\n1o5b9n072m/dd9T69aq6Kqrqq9hdv5uC7AKG5A1pvwAMyh3EwJyBvuR2LAfkDKB/dv8uS05mjqqr\nDqGltYWGlgb2N+9nf/N+6pvqqW+up66pjrqmOuqbOtbrmuqobaplX+M+ahtrqW1qK43+vb2Ne9nb\nsPeA9YaWBgqzCxmQM4ABOQPa/64G5vjl4NzBDMod1F4G53W8HpI3hEG5g9L3x0BtrR/NtnVrR4m8\n3rzZX/i3bvXdoMaN8yXS7/XYY/3Ff9y4PtOTJ5WUCJJlyxbf8+idd3zj1IoV/h9s9K3qiSf6UTtF\nRT1SR9nqWqnZX9NeHbCrfhc1Df5XZM3+GmoaatqXexr2tF+kokttYy2NLY3kZeWR1y/vgGV+Vj45\nmTnk9MshOzP7oPWsjCz6ZfQjK9Mv+2X0a38vMyOTDMsg09qWba8zLAPD2hNPZN3wrx0O59xBy1bX\nSqtrpcW1+GVrS/vrltYWmlubuyxNrU00tTbR2NJIU4tfb2rxrxtaGmhobjhgPbKMXPibW5vJ7ZdL\nbr9c8rLy/LKfPzeRc5Sfld/+XkFWAQXZBRRkFdA/u/8B64U5hX6ZXUhhTiGF2YXkZeWlz6/zw3HO\n3yVHhsBHLysrDy47dvgfUMXFHTM2Rs/eOGqUv+iPHdsrJ23r7ZQIeopz/tdOpP5y5UrfgLV2rf9s\n4kQ/Ycsxx3Ssjx3r/0MMHhy0Mau5tbn9125dUx31zfXtv3wjF8fGlsaDLpqRC237Rbelqf29zhfr\n6NcO13bKOi704JNAdGKIXnaVVCKvI0moc8nMyCQrI4vszGyyMrPIysgiK7PtdUbWQQkuehm56Gdl\nZB09d0tNTX6UWmSkYaREjzyMDFGPHq6+e7e/4Gdn+wEM0RNiDR3q74RHjOiYYTGy3r+/GmnTlBJB\nCLt3+77MGzceuIzMe11ff/AvpsjMdJ1npCsq0q1yX9DUdOAIva5KZGTf3r0HjvSLHpocvd7U1DG3\nSGTekej1QYM6Jq3qXIYMSb8hytJtSgTpqK6uo+502zZf7RQ9I13n2elycw/8zxv5Dz1w4IHzQUSX\nyEimQ835kJ3d+0c4JYNzfjRaU1PHqLbOJTJcOLpE3o8Mle48dDoynLirYdaR4caRocd1db5aJT+/\nY16Prkphof88soxe7zxEubDQ/13rF7qgRND7Oed/9UXfzkcvO88DET03ROcLUPR6Y6Pvdtd5aGz0\nPA6dl5E5HaLnc4ie1yF6LofO64f783UukfcPNVS4peXgEnn/UPNBRA9L7jxMOTOzY56LSInMcRE9\nZLirYcR5eX7ZeTKlSMLtnICjhxxHDz3OytJFW3qMEoF0zTl/Iew8SU70xbLzMnLRjcznEF0iE/VE\nLyPrh4shkig6F+h68iCzrhNRZLtDzQcRSWjRyS2yrjsjOcopEYiI9HHdTQT6iSQi0scFSQRmdpmZ\nrTazdWb2DyFiEBERL+WJwMwygd8AlwEnANeZ2fGpjqMnlZWVhQ4hIYo/nN4cOyj+3irEHcEZwHrn\nXLlzrgl4Gvh8gDh6TG//x6T4w+nNsYPi761CJILRwKao15vb3hMRkQBCJAJ1BxIRSSMp7z5qZmcB\npc65y9pe/yPQ6pz7WdQ2ShYiIt3QK8YRmFk/YA1wEbAVeA+4zjn3UUoDERERAFL+yB/nXLOZfQN4\nFcgE/qAkICISTlqOLBYRkdQJOrI4loFlZvbrts+Xmtn0VMd4OEeK38xKzKzGzBa3lR+FiLMrZvaQ\nme0ws+WH2Sadz/1h40/zcz/WzOab2UozW2Fmdxxiu7Q8/7HEn+bnP9fMFprZkrb4Sw+xXbqe/yPG\nH/f5d84FKfhqofXABCALWAIc32mby4E5betnAn8NFW834y8BZoeO9RDxfwaYDiw/xOdpe+5jjD+d\nz/1I4NS29f74NrPe9G8/lvjT9vy3xZfftuwH/BU4s7ec/xjjj+v8h7wjiGVg2VXAowDOuYXAIDMb\nkdowDynWgXFpOeewc24BsPswm6TzuY8lfkjfc7/dObekbX0f8BEwqtNmaXv+Y4wf0vT8Azjn6tpW\ns/E/5DpPn5u25x9iih/iOP8hE0EsA8u62mZMD8cVq1jid8DZbbeWc8zshJRFl7h0Pvex6BXn3swm\n4O9sFnb6qFec/8PEn9bn38wyzGwJsAOY55x7v9MmaX3+Y4g/rvOf8l5DUWJtpe6c1dKldTuWOD4E\nxjrn6szss8ALwJSeDSup0vXcxyLtz72Z9QeeBb7V9sv6oE06vU6r83+E+NP6/DvnWoFTzWwg8LyZ\nTXPOrey0Wdqe/xjij+v8h7wj2AKMjXo9Fp91D7fNmLb30sER43fO7Y3cwjnn5gJZZlaUuhATks7n\n/ojS/dybWRbwHPC4c+6FLjZJ6/N/pPjT/fxHOOdqgPn4STCjpfX5jzhU/PGe/5CJYBFwrJlNMLNs\n4FpgdqdtZgM3QvuI5Grn3I7UhnlIR4zfzEaY+cdwmdkZ+O66u1Ifarek87k/onQ+921x/QFY5Zz7\n1SE2S9vzH0v8aX7+h5rZoLb1PGAmvp0jWjqf/yPGH+/5D1Y15A4xsMzMbmv7/HfOuTlmdrmZrQdq\ngVtCxdtZLPED1wBfN7NmoA74UrCAOzGzp4DzgaFmtgm4G9/olPbnHo4cP2l87oFzgC8Dy8xscdt7\nPwTGQa84/0eMn/Q+/8XAo+anxM8Anmk7373i2kMM8RPn+deAMhGRPk6PqhQR6eOUCERE+jglAhGR\nPk6JQESkj1MiEBHp45QIRET6OCUCEZE+TolARKSPUyIQ6cTMvmJmT5nZzE7v55lZWWTofqfPss3s\nL22jPUV6FSUCkYPlOeeuc8691un9WcBzrovh+M65RuAN/JxTIr2KEoHIwc40s+O6eP964MXDfO8F\n4IaeCUmk52iuIZEoZnYp8CngXOfcFVHvZwMVzrniw3w3E9jmnBve85GKJI/uCETamNlk4HTn3D3A\n+E4fDwWqD/d951wL0GhmBT0UokiPUCIQ6TALeKptfWOnz+qB3MgLM7vdzBab2YdmFn2XkAPs79kw\nRZIr5KMqRdJNHrCl7UlO66M/cM7tNrNMM8t2zjU6534L/DZ6GzMbAnzSdmcg0msoEYh0uB+4CRgM\n/K8uPp8HfAbfO6grFwAv90xoIj1HjcUiMTKz6cB3nHM3HuLz54B/cM6t7+pzkXSlNgKRGDnnFgPz\nzeyg/zdtD3N/QUlAeiPdEYiI9HG6IxAR6eOUCERE+jglAhGRPk6JQESkj1MiEBHp45QIRET6OCUC\nEZE+7v8D9xyX4J+4fpkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x3f53630>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.1521\n",
      "8.34200492185e-06 9.26558371452e-06 8.90964653776e-06\n",
      "8.45676333368e-06 9.3823368939e-06 9.02458409305e-06\n",
      "8.5712962295e-06 9.49889391734e-06 9.13998574244e-06\n",
      "8.68559293776e-06 9.61524139613e-06 9.25578750397e-06\n",
      "8.79964332369e-06 9.73136668358e-06 9.37192995329e-06\n",
      "8.91343778002e-06 9.8472578486e-06 9.48835789299e-06\n",
      "9.02696721699e-06 9.96290365024e-06 9.60502004587e-06\n",
      "9.1402230518e-06 1.00782935131e-05 9.72186877075e-06\n",
      "9.25319719761e-06 1.01934175037e-05 9.83885979922e-06\n",
      "9.36588205201e-06 1.03082663072e-05 9.95595199182e-06\n",
      "9.47827048509e-06 1.04228312053e-05 1.00731071124e-05\n",
      "9.59035582729e-06 1.05371040541e-05 1.01902896193e-05\n",
      "9.70213185683e-06 1.06510772634e-05 1.0307466472e-05\n",
      "9.81359278706e-06 1.07647437755e-05 1.04246069523e-05\n",
      "9.92473325356e-06 1.08780970454e-05 1.05416824987e-05\n",
      "1.00355483011e-05 1.0991131021e-05 1.06586665534e-05\n",
      "1.01460333706e-05 1.11038401242e-05 1.07755344203e-05\n",
      "1.0256184286e-05 1.12162192315e-05 1.08922631339e-05\n",
      "1.03659972411e-05 1.13282636564e-05 1.10088313379e-05\n",
      "1.04754687864e-05 1.14399691315e-05 1.11252191721e-05\n",
      "1.05845958162e-05 1.15513317904e-05 1.12414081689e-05\n",
      "1.06933755556e-05 1.16623481519e-05 1.13573811564e-05\n",
      "1.08018055477e-05 1.17730151025e-05 1.14731221688e-05\n",
      "1.09098836408e-05 1.18833298806e-05 1.15886163635e-05\n",
      "1.10176079761e-05 1.19932900612e-05 1.17038499445e-05\n",
      "1.11249769752e-05 1.21028935399e-05 1.18188100905e-05\n",
      "1.1231989328e-05 1.2212138519e-05 1.19334848893e-05\n",
      "1.1338643981e-05 1.23210234921e-05 1.2047863276e-05\n",
      "1.14449401253e-05 1.24295472308e-05 1.21619349759e-05\n",
      "1.15508771855e-05 1.2537708771e-05 1.2275690452e-05\n",
      "1.16564548083e-05 1.26455073995e-05 1.23891208552e-05\n",
      "1.17616728519e-05 1.27529426412e-05 1.2502217979e-05\n",
      "1.18665313749e-05 1.28600142473e-05 1.2614974217e-05\n",
      "1.19710306266e-05 1.29667221826e-05 1.2727382523e-05\n",
      "1.20751710363e-05 1.30730666143e-05 1.28394363748e-05\n",
      "1.21789532041e-05 1.3179047901e-05 1.29511297394e-05\n",
      "1.2282377891e-05 1.32846665813e-05 1.30624570414e-05\n",
      "1.23854460099e-05 1.3389923364e-05 1.31734131336e-05\n",
      "1.24881586169e-05 1.34948191181e-05 1.3283993269e-05\n",
      "1.25905169022e-05 1.35993548624e-05 1.33941930752e-05\n",
      "1.26925221825e-05 1.37035317569e-05 1.35040085303e-05\n",
      "1.27941758923e-05 1.38073510936e-05 1.36134359409e-05\n",
      "1.28954795766e-05 1.3910814288e-05 1.37224719207e-05\n",
      "1.29964348834e-05 1.40139228705e-05 1.38311133714e-05\n",
      "1.30970435565e-05 1.41166784789e-05 1.39393574645e-05\n",
      "1.31973074286e-05 1.42190828507e-05 1.40472016242e-05\n",
      "1.32972284148e-05 1.43211378155e-05 1.41546435114e-05\n",
      "1.33968085061e-05 1.44228452885e-05 1.42616810095e-05\n",
      "1.34960497634e-05 1.45242072633e-05 1.43683122098e-05\n",
      "1.35949543119e-05 1.46252258062e-05 1.44745353992e-05\n"
     ]
    }
   ],
   "source": [
    "# As in\n",
    "# Chung, Ting Horng, et al. \"Generalized multiparameter correlation for nonpolar and polar fluid transport properties.\" \n",
    "# Industrial & engineering chemistry research 27.4 (1988): 671-679.\n",
    "\n",
    "from __future__ import print_function\n",
    "import CoolProp, numpy as np, matplotlib.pyplot as plt\n",
    "from math import sqrt, exp\n",
    "%matplotlib inline\n",
    "\n",
    "a0 = [0, 6.32402, 0.12102e-2, 5.28346, 6.62263, 19.74540, -1.89992, 24.27450, 0.79716, -0.23816, 0.68629e-1]\n",
    "a1 = [0, 50.41190, -0.11536e-2, 254.20900, 38.09570, 7.63034, -12.53670, 3.44945, 1.11764, 0.67695e-1, 0.34793]\n",
    "a2 = [0, -51.68010, -0.62571e-2, -168.48100, -8.46414, -14.35440, 4.98529, -11.29130, 0.12348e-1, -0.81630, 0.59256]\n",
    "a3 = [0, 1189.02000, 0.37283e-1, 3898.27000, 31.41780, 31.52670, -18.15070, 69.34660, -4.11661, 4.02528, -0.72663]\n",
    "\n",
    "def check_plot(fluid):\n",
    "\n",
    "    def __G_2(delta):\n",
    "        Y = delta/6.0\n",
    "        G_1 = (1.0-0.5*Y)/(1-Y)**3\n",
    "        G_2 = (A[1]*(1-exp(-A[4]*Y))/Y + A[2]*G_1*exp(A[5]*Y) + A[3]*G_1 )/(A[1]*A[4]+A[2]+A[3])\n",
    "        return G_2\n",
    "\n",
    "    for acentric in [0, 0.489, 0.907]:\n",
    "\n",
    "        A = [0.0]*11\n",
    "        for i in range(1,11):\n",
    "            A[i] = a0[i] + a1[i]*acentric + a2[i]*mu_r**4 + a3[i]*kappa\n",
    "\n",
    "        XXX = np.linspace(0, 3.5)\n",
    "        YYY = [__G_2(delta) for delta in XXX]\n",
    "        plt.plot(XXX,YYY,label=str(acentric))\n",
    "\n",
    "    plt.xlim(0,3.5)\n",
    "    plt.ylim(0,8)\n",
    "    plt.axhline(1.0, dashes = [2, 2])\n",
    "    plt.xlabel(r'$\\delta$ (-)')\n",
    "    plt.ylabel(r'$G_2$ (-)')\n",
    "    plt.show()\n",
    "    \n",
    "check_plot(fluid)\n",
    "\n",
    "fluid = 'n-Propane'\n",
    "Vc_cm3mol = 1/(CoolProp.CoolProp.PropsSI('rhomolar_critical', fluid)/1e6)\n",
    "acentric = CoolProp.CoolProp.PropsSI('acentric', fluid)\n",
    "M_gmol = CoolProp.CoolProp.PropsSI('molemass', fluid)*1000.0\n",
    "Tc = CoolProp.CoolProp.PropsSI('Tcrit', fluid)\n",
    "mu_r = 0\n",
    "kappa = 0\n",
    "print(acentric)\n",
    "    \n",
    "A = [0.0]*11\n",
    "for i in range(1,11):\n",
    "    A[i] = a0[i] + a1[i]*acentric + a2[i]*mu_r**4 + a3[i]*kappa\n",
    "        \n",
    "F_c = 1 - 0.2756*acentric + 0.059035*mu_r**4 + kappa\n",
    "sigma_A = 0.809*Vc_cm3mol**(1.0/3.0)\n",
    "epsilon_over_k = Tc/1.2593\n",
    "\n",
    "def chung(T, rho_molm3):\n",
    "    rho_molcm3 = rho_molm3/1e6\n",
    "    Tstar = T/epsilon_over_k\n",
    "    Omega_2_2 = 1.16145*pow(Tstar, -0.14874) + 0.52487*exp(-0.77320*Tstar) + 2.16178*exp(-2.43787*Tstar);\n",
    "    eta0 = 4.0785e-5*sqrt(M_gmol*T)/(Vc_cm3mol**(2.0/3.0)*Omega_2_2)*F_c\n",
    "    \n",
    "    Y = rho_molcm3*Vc_cm3mol/6.0\n",
    "    G_1 = (1.0-0.5*Y)/(1-Y)**3\n",
    "    G_2 = (A[1]*(1-exp(-A[4]*Y))/Y + A[2]*G_1*exp(A[5]*Y) + A[3]*G_1 )/(A[1]*A[4]+A[2]+A[3])\n",
    "    eta_k = eta0*(1/G_2 + A[6]*Y)\n",
    "    \n",
    "    eta_p = (36.344e-6*sqrt(M_gmol*Tc)/Vc_cm3mol**(2.0/3.0))*A[7]*Y**2*G_2*exp(A[8] + A[9]/Tstar + A[10]/Tstar**2)\n",
    "\n",
    "    # Poise to Pa*s\n",
    "    eta0_Pas = eta0/10.0\n",
    "    eta_k_Pas = eta_k/10.0\n",
    "    eta_p_Pas = eta_p/10.0\n",
    "    return eta0_Pas, eta_k_Pas, eta_p_Pas\n",
    "\n",
    "for T in np.linspace(300, 500):\n",
    "    rho = 1000\n",
    "    eta0, eta_k, eta_p = chung(T, rho)\n",
    "    print(eta0, eta_k + eta_p, CoolProp.CoolProp.PropsSI(\"V\",\"T\",T,\"Dmolar\",rho,fluid))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n"
   ]
  }
 ],
 "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
