CoolProp/doc/notebooks/braindump/iPRSV.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The iPSRV from Teus van der Stelt et al.\n",
      "$$\\kappa = \\kappa_0+\\kappa_1\\left\\lbrace \\sqrt{[A-D(T_r+B)]^2+E}+A-D(T_r+B) \\right\\rbrace\\sqrt{T_r+C}$$\n",
      "$$\\kappa_0 = 0.378893 + 1.4897153\\omega-0.17131848\\omega^2+0.0196554\\omega^3$$\n",
      "$$ p = \\frac{RT}{v-b}-\\frac{a}{v^2+2 b v-b^2}$$\n",
      "with\n",
      "A = 1.1\n",
      "B = 0.25\n",
      "C = 0.2\n",
      "D = 1.2\n",
      "E = 0.01\n",
      "\n",
      "For ethanol, the value $\\kappa_1$ is -0.03374\n",
      "\n",
      "For water, the value $\\kappa_1$ is -0.06635"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from CoolProp.CoolProp import Props\n",
      "from math import sqrt\n",
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "import scipy.optimize\n",
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 64
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "kappa_dict = dict(Ethanol = -0.03374,\n",
      "                  Water = -0.06635,\n",
      "                  Propane = 0.03161)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 65
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Pure Fluids"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "T = 300 #[K]\n",
      "p = 3000 #[kPa]\n",
      "Fluid = 'Propane'\n",
      "rhoc = Props(Fluid,'rhocrit') #[kg/m^3]\n",
      "omega = Props(Fluid,'accentric') #[-]\n",
      "R = 8.314472/(Props(Fluid,'molemass'))\n",
      "Tc = Props(Fluid,'Tcrit')\n",
      "pc = Props(Fluid,'pcrit')\n",
      "T_r = T/Tc\n",
      "kappa_1 = kappa_dict[Fluid]\n",
      "rho_EOS = Props('D','T',T,'P',p,Fluid)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 66
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Ts = np.linspace(Props(Fluid,'Ttriple'),Tc,300)\n",
      "ps = Props('P','T',Ts,'Q',1,Fluid)\n",
      "plt.plot(Ts,ps)\n",
      "plt.xlabel('T [K]')\n",
      "plt.ylabel('p [kPa]')\n",
      "plt.yscale('log')\n",
      "plt.plot(T,p,'o')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 67,
       "text": [
        "[<matplotlib.lines.Line2D at 0x871c370>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAESCAYAAADqoDJEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH7BJREFUeJzt3XmYVOWZ9/Fvg6zi4BgmCoJToCBLWAQFFcFWUFmMxGhU\nEMfo4GhEcDJxQ983QoKIxoxENEheQAER0IgoINJELNPKLrIvArGRRQGXUZBFoev9466aqm6qu6u7\nluecqt/nus5VXadruTnd9F3PeoOIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIh4TnXXAVTgZGAx\nsAfY6jgWEREBqrkOoAIPADNdByEiIlGZThyTgL3AulLnewGbsVbFg+FzVwAbgf0Zi05ERDynG3Ae\nJRNHdWAbEABqAKuBVsBI4GlgATAbyMtkoCIi4h0BSiaOi4C3Y+4/FD4ibgX6pD8sERFJxEmuAwDO\nBHbG3N8FdIm5Pzmz4YiISHm8kDhCVX1io0aNQnv27EllLCIiuWAP9qG9Srwwq2o30CTmfhOs1VGh\nPXv2EAqFfHs8+uijzmNQ/O7jUPz+O/wceygUAmiUzB9tLySOlUBzbOyjJnAj8KbLgEREpGyZThzT\nsQV9LbBxjduAY8A92Oypjdi6jU0ZjktERBKU6TGO/mWcnx8+ckp+fr7rEJKi+N1S/O74OfZU8Pva\niFC4v05ERBKUl5cHSfz998IYh4hIzjl4EFasgFdecR1J5XlhOq6I5JB5C+fxzMvPcDR0lFp5tRg6\nYCh9r+jrOqy0OXwYNm2CDRvsWL/ebvfuhZYtoW1buP56qOajj/FKHCKSMfMWzuPe5+5l+3nb//fc\n9ufsa78nj6NHYcuWaGKIHLt2QfPm0KYN/OQnMGiQfd2sGVT3+v7kZdAYh4hkzFW3XUVBoODE8zuu\n4u1Jb8d5hvf88ANs3RpNEJHbHTugaVNLDm3aRI9zzoEaNVxHXVKyYxxqcYhIxhwNHY17/kjxkQxH\nkpj9+2HtWjvWrLHbzZuhceNogvjFL2DECGjRAmrWdB1xZihxiEjG1MqrFfd87Wq1MxxJST/8YN1M\nkeQQuT10CNq3h3bt4JJL4O67LVmcfLLTcJ1T4hCRjBk6YCjbn9teYozj7FVnM+SeIRmL4euv4aOP\nYPXqaILYsgXOOssSRPv2MHiw3TZpAnl+79BPA79fEo1xiPjMvIXzGDt9LEeKj1C7Wm2G9B+StoHx\n/fth1So7PvzQbr/4wpLCeedFWxNt2kDdumkJwZOSHePwcuJoCdwL/AjbjmRinMcocYgIAJ99VjJB\nrFoF334LHTva0amT3TZv7q+pr+mQzYkjohowA7ghzveUOERy0JdfwvLldqxYYcni++8tOUQSRMeO\nNuVVXU0n8lvimAT0BfYBbWPO9wLGYGVkJwBPhM//FLgb+H/ArDivp8QhkuWOHLHxiOXLYdkyO/bv\nh/PPh86d4YILLFmcdZaSRKL8lji6AQeBKUQTR3VgC9ATq82xAtsMMXaH3DeAfnFeT4lDJIsUF9sa\nidgksXEjnHsudOliiaJLF7vv18VzXuC3dRyFWN2NWJ2BbUBR+P4MLEn8GPg5UBt4NzPhiUgmHT5s\nSeKDD+xYsgTq17fk0KUL9O9vg9i5NHDtB16YjltWzfH3woeIZIl9+yxBvP++3a5fbwvpLrnEtuKY\nNAlOP911lFIRLySOpPqa8vPzCQQCBAIB8vPzc36ffBGvCIVsfcT770cTxZdfwkUXQdeu8OSTNj5R\np47rSLNfMBgkGAxSVFREUVFR0q/nYigpAMwhOsZxITAcGyAHGAYUEx0gL4/GOEQ8IhSCjz+GYBDe\nfddua9eG7t0tUXTtCq1bayqsF/htjCOe2Jrje7Ca42VVChQRjwiFYNu2komiRg247DLo1QtGj4ZA\nwHGQkhaZbnFMBy7FFvXtA34LvAD0JjoddyLweIKvpxaHSAZ9+iksXBhNFHl5lijy8+02ENCUWD/w\n23TcVFPiEEmjgwctQRQUWML44gvo2RN69LBEoQV2/qTEocQhkjLHj9tWHZFE8eGHNoB95ZV2dOig\nMYpsoMShxCGSlC++gPnzYd48SxZnnBFNFN27awvxbKTEocQhUimhkK2fmDvXjvXr4fLLoW9fG9Ru\n3Nh1hJJuShxKHCIVOnzYxioiyaJaNbj6ajsuvdSmzUruyIbpuCKSBl99BXPmwKxZNguqQwdLFPPn\nQ6tWGtSWqvP7r45aHCIxPv8cZs+2ZLF0qc1++vnPrRvqtNNcRydeoa4qJQ7JcUVF8PrrlizWr4fe\nveG662y8QgPbEo8ShxKH5KBdu2DmTJgxwxLHNddYsujRA2rVch2deJ0ShxKH5IgvvoC//hWmT4d1\n6+Daa23b8fx8OEmjlVIJShxKHJLFDhyAN96wZPH++9YN1b+/dUOpZSFVlc2Jox9WZvafsP2rFsZ5\njBKHZJ3iYli0CF580abOXnIJDBhg3VH16rmOTrJBNieOiFOBp4BBcb6nxCFZY9s2mDzZjgYN4Je/\ntNbFv/yL68gk2ySbODK968wkYC+wrtT5XsBmYCvwYKnv/R/g2fSHJpJ5Bw5Y1bvu3eHii21TwTlz\nbL+ooUOVNMSbMt3i6AYcBKYQLeRUHdgC9AR2AyuwehybgdFAAfBOGa+nFof40sqV8PzzNtidnw+3\n3WbjFzVruo5McoHfVo4XYgWbYnUGtgFF4fszsPGNnkAPbIzjHGB8RiIUSZPvvrNB7ueft1Xdd95p\npVVVY1v8xguT+M4Edsbc3wV0AYYAYyt6smqOi9etWwfjx1vS6N4dHnsMrrhC25NL5qS65rgXEkdS\nfU3BYDBFYYikzvHj8OabMGYMbN8OgwbBmjXaeVbcKP2hOi/Jjcq8kDh2A01i7jfBWh0ivvPNNzbY\nPXas1bX4z/+0vaK0QE+yiRd+nVcCzbGxjz3AjdjguIhvbN9uyWLKFLjqKuuW6tLFdVQi6ZHpXtbp\nwGKgBTaucRtwDLgHWABsBGYCmzIcl0iVLFtmW39ceCHUqQNr1yppSPbzwwLA8mg6rmRcKATvvAOP\nP26L9u6/36bTaida8Qu/TccV8a3iYts3atQoW6j30EO2FUiNGq4jE8ksJQ6RChw7Bi+/DKNHW6vi\n4YehXz9Np5XcpcQhUobjx63exYgR0LAh/OlP0LOnSq6KKHGIlFJcbNX0Hn0U6teHcePg8suVMEQi\nlDhEwkIh22Dwt7+1cYunnrK6F0oYIiUpcYgAixfDfffZflK/+53VvlDCEIlPiUNy2tatNjtqxQrb\nQ+rmmzXoLVIR/ReRnLR/PwwZAhddBBdcYLvU3nKLkoZIIvTfRHLK0aPwxBPQqpV1RW3aZC2OOnVc\nRybiH17uqmoKPALUB37hOBbJAm+/bVX1WraEJUugeXPXEYn4kx+G/16l7MShLUekQp98Ar/+NWzY\nYGsx+vRxHZGIW7lQc1ykSg4fhuHDbQyjc2dYv15JQyQVMp04XsCSRKzqwLPh862xLdVbZTguyTJv\nvw1t2lgrY9Uq2yakVi3XUYlkh0wnjkLg61LnYmuO/0C05vhpwPNAB9QKkQTt3w8DB8Ldd1tt71df\nhbPOch2VSHbxwuB4WTXHvwLuqujJqjkuYKu+p02D3/zGEse6ddrmXCRCNcdLUc1x2bED7roL9uyB\nuXNtTENEolJdc9wL6zhUc1yqJBSy+t7nnw/dusHKlUoaIpnghRaHao5Lpe3dC3fcAZ9+CosWQdu2\nriMSyR2qOS6+M2sWtG9vyWL5ciUNkUzzwwLA8mgBYA759lu45x5b9T1liu0zJSKV57cFgCJV8uGH\n0KmT7Sm1erWShohLShziaaEQjBljBZUeewzGj9c0WxHXvDA4LhLXl1/C7bfbNNtly6BZM9cRiQio\nxSEetXgxnHcenH02fPCBkoaIl6jFIZ4SCsG4cbY54cSJ8NOfuo5IREpT4hDPOHLE9phavtxaHOec\n4zoiEYlHXVXiCTt3QvfucOAALF2qpCHiZUoc4tx771m9jOuvh1degXr1XEckIuVRV5U4NWECPPII\nTJ0KV17pOhoRSYSXE8fJwJ+Bo0AQeNlpNJJSxcUwbBi89hr8/e9w7rmuIxKRRHl5y5FbsJoc87Di\nTjfFeYy2HPGhQ4fg3/7NNip8/XVo0MB1RCK5xW9bjlSm5nhsgafjGYlO0u7zz+Gyy6B2bfjb35Q0\nRPzIyzXHdxGt06FB/CywfTt07Wrbh0ydqhrgIn7l5Zrjs4DrsHGONzMXoqTDRx9ZsaUHHoARIyDJ\nAmQi4pAXBsfLqjl+CLjdSUSSUsEg3HCDrQi/7jrX0YhIsryQOJIa3c7PzycQCBAIBE6oqyvuzZpl\n9cBnzrSxDRHJvGAwSDAYpKioiKKioqRfz0WHQQCYA0Tqtl0IDCc69jEMKAaeSOC1NKvKw6ZMgYce\ngrlzoWNH19GISITfZlXFE1tzvCZWc1xjGj43cSI8/LDVA1fSEMkuqjkuKff88zYA/u670LKl62hE\nJNX8PrdFXVUe88wz8N//bS0N1dAQ8aZku6q8MDguWeLpp+HZZ23Twn/9V9fRiEi6KHFISowbB2PH\nWtJo0qTix4uIfylxSNKmTIFRo5Q0RHKFxjgkKa+9BkOG2JiGBsJF/EFjHOLMW29ZqdcFC5Q0RHKJ\nWhxSJYWFtn3InDnQpYvraESkMpJtcShxSKVt3Gjbh0ybBj17uo5GRCorG1aOi4/s3g29e8Mf/6ik\nIZKrlDgkYd98Y0lj8GAYONB1NCLiSnlNldJV+uLZD1yeoliqQl1VGXL0qCWNNm1sdbjqaYj4Vzpn\nVVUHelfw4uncjLAp8AhQH/hFGt9HKhAKwb//O/zzP8OYMUoaIrmuvMTxH8COCp4/OIWxlPYJMAh4\nNY3vIQkYPRo+/tgW+FWv7joaEXGtvDGO9xN4fmECj5kE7OXErq9ewGZgK/BgAq8jDrzxBvz5zzB7\nNtSp4zoaEfGCRAbHWwB/xbY6/yR8/KMS7/EC0SJNEdWBZ8PnWwP9gVbALcDTQKNKvL6kybp1cMcd\nVsWvkX4iIhKWSOJ4AXge+AHIByYD0yrxHoXA16XOdQa2AUXh150B9AOmAr8G9gCnhd+3A2qRZNz+\n/dCvn41pXHCB62hExEsS2XKkDvA3bJB8B1bmdRXwf5N43zOxQk4Ru4DS64+/Au6q6IVUczz1vv8e\nrr8ebroJBgxwHY2IJCvVNccTSRxHsK6lbVilvj3AyUm+b8rm0AaDwVS9lITdfz/Urw8jR7qORERS\nofSH6rwkp0YmkjjuBeoCQ4HfA/8E3JrUu8JuIHYD7iZYq0Mce+UVmDsXPvwQqml5qIjEUV7iOB14\nGDgHWAs8DvwyRe+7EmgOBLAWzI3YALk4tGWLrQpfsABOPdV1NCLiVeV9ppwCHATGAqcAz1TxPaYD\ni7HZWTuB24BjWLfXAmAjMBObtSWOfPed7XY7ahR07Og6GhHxsvI6utYA7WPufwScl95wKk1bjqRA\nKAS33morwl98USvDRbJdOrccycOmxEa+rh5zH2zWk2SBCRNg1SpYtkxJQ0QqVt6fiSLKn/3UNLWh\nVIlaHEnasAHy860wk6r4ieSGdLY4AuW9b1XfULzjyBFbpzF6tJKGiCQukQmXvyt1vzrwUhpikQx7\n+GFo3hxuv911JCLiJ4kkjrOAYeGvawGzsI0JxccKCuDVV+Evf9G4hohUTiJ/Mqphe1OtxYo2vYVt\nROgFGuOogv37oUMHmDoVLndZhktEnEh2jKO8J3YiOjheAxiPrceYED63qqpvmkJKHJUUCtk+VM2a\nwR/+4DoaEXEhnYkjSMlZVXml7l9W1TdNISWOSpoxA37/e9tSpHZt19GIiAvpTBx+oMRRCXv3Qrt2\nMGcOdO7sOhoRcSXZxFHe4PjVCTw/kcckox/wF6xexxVpfq+sFgrBr35lM6iUNEQkGeVlnM3AAE7s\noop97otA29SHdYJTgaewGuSx1OJI0PTptk36qlVQq5braETEpUyOccTzFXBdAu8zCegL7KNkoukF\njMHWhkwAnijj+U9ha0dWlzqvxJGAffugbVvbLl3V/ETEL2Mc3bCddqcQTRzVgS1AT6w+xwpsa/Xz\ngY7AH4DPgNFAAfBOnNdV4kjAwIHQsKFmUYmISeeWI6lUyIlbmMTWHYdo3fHRWO1xsOJRPbDiUedg\nU4KlEhYuhA8+gPXrXUciItkiU4kjnkTqjj9DBXVAVHO8bIcP24D4c8/ByckW+xUR33JRczxdUtLH\npJrjZRs50ooy9enjOhIRcclFzfE6wN3AJdgf+0JgHHAkqXdW3fG02rDB9qFas8Z1JCKSbRLZ5HAK\n0BrrMnoWaEN0DCIZsXXHa2J1x99MwevmvMiajREjoFEj19GISLZJpMXRBkscEYuwOuGVMR24FPgR\nNq7xW+AFonXHqwMTUd3xlJgxAw4ehDvvdB2JiGSjRDq6XgKeA5aE718IDAZuSVdQlaDpuKUcPAit\nWlny6NrVdTQi4kWZWMexGWiBtRRCWH2OLcCx8P12VX3zFFDiKOWRR2DHDnhJpbZEpAyZSByBCr5f\nVNU3TwEljhjbt0OXLjYgfuaZrqMREa/yy8rxdFHiiPGzn1niGDas4seKSO7yy8pxSbOCAlsdPmOG\n60hEJNslMh1XPO74cbj/fnjySRVnEpH0U+LIAtOm2ZYi117rOhIRyQUa4/C5I0fg3HPh5Zc1/VZE\nEpPOCoDiA2PH2n5UShoikilqcfjYV19Za6OwEFq2dB2NiPhFNk/HbQnci21TsgDbkqS0nE4c990H\nBw7AeFUpEZFKyObEEVENK/J0Q5zv5Wzi2LkTOnSwKbgNG7qORkT8xA9jHJOAvcC6Uud7YduZbAUe\nLOO5PwXmYYlDYowaBYMGKWmISOZlosVR1Xrje2Je4w2srGxpOdni2LHDBsS3bIEGDVxHIyJ+44eV\n41WtN34p8HOgNvBuuoP0k8ceg7vuUtIQETdcbTmSSL3x98KHxPjHP2DWLPj4Y9eRiEiucpU4Uta/\nlJ+fTyAQIBAInFBXNxuNHAmDB8Npp7mORET8IhgMEgwGKSoqoqioKOnXy9SsqgAwh+gYx4XAcGyA\nHGAYUAw8UcnXzakxjq1b4aKLYNs2OPVU19GIiF/5YVZVPKo3XgUjR8LQoUoaIuJWJrqqVG88BXbs\ngLlzrViTiIhLflgAWJ6c6aoaMgTq1oUnKtuZJyJSSi6sHC9PTiSOfftsL6oNG7TgT0SS59cxDqmE\nZ56BG25Q0hARb1CLw+O+/RaaNYNly+Dss11HIyLZQC2OLDd+PFx5pZKGiHiHWhweduSItTbmz4f2\n7V1HIyLZQi2OLDZtmm2drqQhIl7iassRqUAoBGPGwNNPu45ERKQktTg8atEiSx49eriORESkJCUO\nj/rTn+DeeyHP76NQIpJ1vJ44TsaKPPV1HUgmbd0KS5fCwIGuIxEROZHXE8cDwEzXQWTa2LFWFrZO\nHdeRiIicyMs1x68ANgL70xqdx3zzDbz0Etx9t+tIRETiy8SsqheAsVjN8YjqwLOUrDn+JiVrjl+K\ndVW1Bg4Db5HCAlBeNWkS9OoFjRu7jkREJD5XhZwuAh4lWsjpofDt6DjPvRVrdbwV53tZtQDw+HFo\n3hymT4cupQvpioikSLILAL1cczxicvrD8YaCAisJq6QhIl6mmuMeMn483HWX6yhEJNuo5nhJWdNV\ntWsXtGsHn34K9eq5jkZEsplf96pSzfFSJk6Em25S0hAR78tE4pgOLAZaYOMatwHHiNYc34it1cjZ\nmuPHjsGECXDnna4jERGpmN83tMiKrqo5c2DUKFiyxHUkIpIL/NpVJTHGj1drQ0T8Qy0Ox3bsgI4d\nYedOqFvXdTQikgvU4vC5iRPh5puVNETEP9TicKi4GJo2tTGOdu1cRyMiuUItDh8LBm2luJKGiPiJ\nEodDkyfDrbe6jkJEpHLUVeXIwYPQpAls2QI//rHraEQkl6iryqdeew26dVPSEBH/UeJwRN1UIuJX\n6qpyYMcO6NQJdu+GWrVcRyMiuSabu6rygUJgHFYNMGtMmQI33qikISL+5KoeRyKKgQNALazQU1YI\nhSxxTJvmOhIRkarJRItjErAXWFfqfC9gM7AVeDDO8wqBPlhZ2RHpDDCTli6Fk06CCy5wHYmISNVk\nInG8QLRgU0R14Nnw+dZAf6AVcAvwNNCIaJXA/8FaHVlhxgzo3x/y/D66JCI5KxNdVYVYwaZYnYFt\nQFH4/gygHzAamBo+dy1wFXAqMDbdQWbC8ePwyiu2YlxExK9cjXGciRV1itgFdCn1mNfDR7n8VHO8\nsBDOOAPOPdd1JCKSS1Jdc9xV4kjZHNqgjz6+z5hh5WFFRDKp9IfqvCT7yl0ljt1Ak5j7TciimVPx\n/PCDrRZfvtx1JCIiyXG1jmMl0Bwb+6gJ3Ai86SiWjFi0CM4+27ZRFxHxs0wkjunAYqAFNq5xG3AM\nuAdYAGwEZgKbMhCLMzNn2qI/ERG/8/ukUF9sOXL0KDRqBGvWQOPGrqMRkVyXzVuOZI2CAmjTRklD\nRLKDEkcGqJtKRLKJuqrS7OhRW7uxaZPdioi4pq4qj3vnHfjJT5Q0RCR7KHGk2euvw7XXuo5CRCR1\n1FWVRseP22yqJUugWTPX0YiIGHVVedjixdZFpaQhItlEiSON1E0lItnIyxUAfS0UssQxe7brSERE\nUsvLiSMPGAmcgu1tNcVtOJWzdq0Va2rXznUkIiKp5eWuqp9hdTu+x4c750a6qVTpT0SyjZdrjrcA\nPgDuA36VzgDTYc4cuOYa11GIiKSel2uO78LqjQMUZyDOlNmzBz75BC6+2HUkIiKp5+Wa47OwWuPd\ngGCaY0ypt96Cq66CGjVcRyIiknperjl+GBiUsYhSaO5cuP5611GIiKSH72uO5+fnEwgECAQCJ9TV\ndeHIEav2N2GC0zBERP5XMBgkGAxSVFREUVFR0q/n+5rjwWAwFfGkzHvv2RTcBg1cRyIiYkp/qM5L\ncrqnao6n2Ny50Lev6yhERNJHNcdTKBSyxHH11a4jERFJH78vT/PU7rgbNkCfPlBUpIV/IuJd2h3X\nQ+bPh969lTREJLspcaRQQYGt3xARyWZ+/2zsma6qQ4fg9NNh1y6oX991NCIiZVNXlUcUFkKHDkoa\nIpL9lDhSZMECuPJK11GIiKSfEkeKaHxDRHKFxjhSYNcuaN8e9u2D6tVdRyMiUj6NcXhAQQH07Kmk\nISK5QYkjBdRNJSK5xMtdVZcAN2MbMbYGusZ5jPOuquPHbRru6tXQuLHTUEREEpJsV5Wr3XET8X74\n6AcsdxxLmQ4dgt/8RklDRHKHl2uORwwAXk5PaMk75RQYNqxqz/XalvCVpfjdUvzu+Dn2VPByzXGA\ns4BvgO8yEGfG+f2XT/G7pfjd8XPsqeDlmuMAt2MtFhER8Qgv1xwHGJ6RaEREJGGZmlUVAOYAbcP3\nr8O6qe4I3x+IJY4hlXzd3US7tUREJDF7sA/wVeL3muNV/oeLiIi3BSg5q+okYDvRmuOrscFxERER\npmPNoqNEa44D9Aa2YIPkVZzQKiIiUr4iYC3wEdFFgacBC4GPgQLgVCeRxRdvDUt58Q7D1rVsBlxv\n0h4v9uFYl+JH4aN3zPe8FDtY9+e7wAZgPTA0fN4v17+s+Ifjj59BbWAZ1puwEXg8fN4v17+s+Ifj\nj+sPtuzhI2x8Gfxz7VPuE+wfH+tJ4IHw1w9iU3q9ohtwHiX/+JYVb2vsl7QG1oW3Dbd7icWL/VHg\nv+I81muxA5wBdAh/XQ9r3bbCP9e/rPj99DOoG749CViKbSPkl+sP8eP30/X/L2Aa8Gb4fsquvet/\nWFWUngl2DTA5/PVk4GeZDadchcDXpc6VFW8/rFvvB6xltQ1b7+JKvNgh/kw8r8UO8Dn2nwHgILAJ\nm0zhl+tfVvzgn5/BofBtTezT79f45/pD/PjBH9e/MdAHmEA03pRde78ljhDwN2Al0am8p2NdKoRv\nT3cQV2WUFW8jSs4s24U3Z40NAdYAE4k2db0eewBrPS3Dn9c/gMW/NHzfLz+Daljy20u0281P1z9e\n/OCP6/80cD9QHHMuZdfeb4mjK/YfqDcwGOtOiRUKH35RUbxe+7eMA5piXSifAX8s57Feib0e8Bpw\nL3Cg1Pf8cP3rAX/F4j+Iv34GxVicjYHuwGWlvu/16186/nz8cf2vBvZh4xtlrdVL6tr7LXF8Fr7d\nD7yONaf2Yv3BAA2xC+ZlZcVbem1L4/A5L9lH9BduAtHmrFdjr4EljanA7PA5P13/SPwvEY3fbz8D\nsP3m5gGd8Nf1j4jEfz7+uP4XY91Sn2BdUJdj/wf8eO2TVhc4Jfz1ycAH2Oj/k0R3130Ibw2Ow4lr\nWMqKNzJAVRP7RLMd9/VSApSMvWHM178mumuxF2PPA6ZgTfZYfrn+ZcXvl59BA6LdOHWAvwM98M/1\nLyv+M2Ie4+XrH3Ep0VlVfrn2KdUU+8etxqYnRtZ+nIaNe3hxOm5kDcv3RNewlBfvw9jA1GbAdU3B\n0rHfjv0hW4v1786m5HiSl2IHmwFTjP2+RKZO9sI/1z9e/L3xz8+gLbAKi38t1t8O/rn+ZcXvl+sf\ncSnRWVV+ufYiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIpX2I6KL8T4jWpdhFbY9SMTx8LnIKuMi\nouUAOgH/ANoDN2I1EOYgIiJZr6y6DHDiZoqROjLtsKRxfsz3YreEEPEFv21yKOIlldnPpw22MedA\nrCxAVV5DxBNOch2ASA7Iw/Y1uhlY7DgWkaSpxSGSfiGs1vMd6P+cZAH9Eotkxj3h2z87jUIkBZQ4\nRDKjGBgAtARGOI5FJClKHCJVl2hp0MjjjmKV2a4BflXJ1xARkSxWejpuefLRdFzxGbU4RFLvW2wB\nYMMKHncj8BzwVdojEhEREREREREREREREREREREREZE0+f8uypooNMv09gAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x83b4770>"
       ]
      }
     ],
     "prompt_number": 67
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "A = 1.1\n",
      "B = 0.25\n",
      "C = 0.2\n",
      "D = 1.2\n",
      "E = 0.01\n",
      "kappa_0 = 0.378893 + 1.4897153*omega-0.17131848*omega**2+0.0196554*omega**3\n",
      "kappa = kappa_0+kappa_1*(sqrt((A-D*(T_r+B))**2+E)+A-D*(T_r+B) )*sqrt(T_r+C)\n",
      "alpha = (1+kappa*(1-sqrt(T_r)))**2\n",
      "a = 0.457235*R**2*Tc**2/pc*alpha\n",
      "b = 0.077796*R*Tc/pc"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 68
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "A = a*p/(R*R*T*T)\n",
      "B = b*p/(R*T)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 69
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Z = np.linspace(0,1)\n",
      "def f(Z):\n",
      "    return 1*Z**3+(-1+B)*Z**2+ (A-3*B*B-2*B)*Z -A*B+B**2+B**3\n",
      "plt.plot(Z,f(Z))\n",
      "plt.axhline(0,linestyle = '--')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 70,
       "text": [
        "<matplotlib.lines.Line2D at 0x8909b50>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEACAYAAAC08h1NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGZFJREFUeJzt3XuUVNWdt/GnA2jUqGiYpYkyoqByCcRgQPASW0WFSSJG\nTRwwXqJGvOD7rjVvVhzUhJ5M3nhZmqCD8YJ4d4LEAEFfFVCoaFTuF0FohB6IgEoMiCCgAlXvH7va\nbttuqO5TXaeqzvNZ66yu6tp9zvaAX3b/zj77gCRJkiRJkiRJkiRJkiRJkkrcQKAaWAHc0MjnlcCH\nwILsdnPBeiZJyrs2wEqgE9AOWAh0a9CmEphc0F5Jkpr0pYg/35cQ/KuBHcA4YHAj7SoiHkeSlCdR\ng/8wYE2992uz36svA5wILAKeA7pHPKYkKYK2EX8+k0Ob+UBHYBswCJgEHBPxuJKkFooa/OsIoV6r\nI2HUX9+Weq+fB34PHAxsrN+oc+fOmZqamojdkaREqQG6NPeHopZ65gJHEy7u7gVcyBcv5B5CXY2/\nb/b1xgZtqKmpIZPJJH4bOXJk7H0ols1z4bnwXOx+Azq3JLijjvh3AsOBKYQZPmOBZcCw7Of3AxcA\n12TbbgP+NeIxJUkRRA1+COWb5xt87/56r+/JbpKkIhC11KM8q6ysjLsLRcNzUcdzUcdzEV0xza/P\nZGtWkqQcVFRUQAty3BG/JCWMwS9JCWPwS1LCGPySlDAGvyQljMEvSQlj8EtSwhj8kpQwBr8kJYzB\nL0kJY/BLUsIY/JKUMAa/JCWMwS9JCWPwS1LCGPySlDAGvySVoA8+aPnPGvySVIJ+8YuW/6zBL0kl\nZv58ePrplv+8wS9JJSSdhuuug9/8puX7MPglqYQ88kj4etllLd9Hs5/O3ooymUwm7j5IUtHauBG6\nd4fnnoPevaGiogJakOMGvySViGuugTZtYPTo8L6lwd82v92SJLWGuXNh4kRYtiz6vqzxS1KRS6fh\n2mvh1lvhoIOi78/gl6QiN3YstGsHl1ySn/3lI/gHAtXACuCG3bTrA+wEzsvDMSUpETZsgJtvhnvu\ngS/laage9eJuG2A5MABYB8wBhgANq1BtgGnANuBh4E+N7MuLu5LUwLBhsPfecPfdX/wsrou7fYGV\nwOrs+3HAYL4Y/NcDTxNG/ZKkHLz2GjzzDCxdmt/9Rv3F4TBgTb33a7Pfa9hmMHBv9r3Deknag08/\nhZ/+FEaNgvbt87vvqMGfS4iPAv4927aC4rp3QJKK0m23wVFHwQ9/mP99Ry31rAM61nvfkTDqr+94\nQgkIoAMwCNgBTG64s6qqqs9eV1ZWUllZGbF7klR6qqtDTX/+fKioN1ROpVKkUqnI+486+m5LuLh7\nBvAOMJvGL+7Wehh4BpjQyGde3JWUeOk0VFaGkf711+++bVwXd3cCw4EphJk7YwmhPyz7+f0R9y9J\nifLgg6G+f+21rXeMYqq3O+KXlGjvvgu9esH06dCz557bu0ibJJW4Cy6Arl3h17/Orb2LtElSCfvz\nn2HxYnjiidY/liN+SYrZ5s3Qo0cI/VNPzf3nLPVIUom67rpwQXfMmOb9nKUeSSpBM2bUlXkKxWWZ\nJSkmW7bA5ZfDAw/kZ539XFnqkaSYXH017NgR1ttvCUs9klRCpk6F55+HN94o/LENfkkqsE2b4Mor\nw0j/wAMLf3xLPZJUYD/5CXz5y3DvvXtuuzuWeiSpBDz7LPzlL/GUeGo54pekAtm4MazB8+STYQXO\nqLyBS5KK3EUXQYcOcNdd+dmfpR5JKmITJsDs2bBoUdw9ccQvSa3unXegd+8Q/ieemL/9tnTE7527\nktSK0mm49NJws1Y+Qz8Kg1+SWtHvfgfbtsHNN8fdkzqWeiSplSxYAGedFWr7Rx6Z//1b6pGkIrJ1\nKwwZEmbwtEboR+GIX5JawbBhocTz+OOtdwync0pSkZg4EaZNg4UL4+5J4xzxS1IerVsXpm5OmgT9\n+7fusazxS1LM0mm45BIYPrz1Qz8Kg1+S8uT228Ozc2+8Me6e7J41fknKg5dfhlGjYM4caNMm7t7s\nniN+SYpo/XoYOhQeeQQ6doy7N3vmxV1JimDXrnCTVv/+8OtfF/bYXtyVpBj86leQycB//EfcPcmd\nNX5JaqGpU8Nzc+fOLf66fn35GPEPBKqBFcANjXw+GFgELADmAafn4ZiSFKu1a8PUzSefhEMPjbs3\nzRO1xt8GWA4MANYBc4AhwLJ6bfYDtmZf9wQmAl0a2Zc1fkklYceO8OjE730PRoyIrx9x1fj7AiuB\n1cAOYBxhhF/f1nqvvwL8I+IxJSlWI0ZA+/ZwQ2M1jhIQtcZ/GLCm3vu1wAmNtDsXuAX4GnBWxGNK\nUmwmTICnn4Z58+BLJTo9Jmrw51qbmZTdTgEeB45trFFVVdVnrysrK6nMx2PoJSlPFi8OT9J67jn4\n6lcLf/xUKkUqlYq8n6g1/n5AFeECL8AIIA3ctpufqSGUiDY0+L41fklFa+NG6NMnTN+86KK4exPE\nVeOfCxwNdAL2Ai4EJjdo07lex3pnvzYMfUkqWjt3woUXwnnnFU/oRxG11LMTGA5MIczwGUuY0TMs\n+/n9wPnAJYSLvx8B/xrxmJJUUD//eajn33pr3D3JD5dskKTdeOwx+M//DM/NPeiguHvzeS0t9Rj8\nktSEOXPgu9+FGTOgR4+4e/NFrtUjSXn03nuhpj9mTHGGfhQGvyQ18MkncP758NOfwuCGt6SWAUs9\nklRPJgMXXwwffwzjxxf3TVotLfW4Oqck1TNyJKxcGer6xRz6URj8kpT1yCPwxBMwcybss0/cvWk9\nlnokCZg+HYYMgVQKunWLuze5cVaPJLXQ0qUh9J96qnRCPwqDX1KivfdemKt/xx1hjf0kMPglJda2\nbXDOOXDZZWEmT1JY45eUSLt2wQUXwAEHhIu6FcWUhjlyOqck5SiTgWuugc2bQ12/FEM/CoNfUuLc\neCMsXAgvvQR77RV3bwrP4JeUKHfcAZMmwSuvwP77x92beBj8khLj4Ydh9OgQ+h06xN2b+BRTZcuL\nu5JazaRJoa6fSsGxjT71u/R4cVeSmjBjBlx1FTz/fPmEfhTO45dU1ubNC8/LHT8ejj8+7t4UB4Nf\nUtlavBi+9z144IHk3JWbC4NfUll680046ywYNQrOPTfu3hQXg19S2Vm6FM48E+68M5R59HkGv6Sy\nUl0dQv/222Ho0Lh7U5wMfkll4623YMAA+M1v4Mc/jrs3xcvgl1QWVqyAM86AX/0KLr007t4UN4Nf\nUsmrqQmhP3IkXH553L0pfga/pJK2bFmYqnnTTXDllXH3pjR4566kkjV/fnh61q23Wt5pDoNfUkn6\n61/hvPPgvvvCV+UuH6WegUA1sAK4oZHPLwIWAW8ArwK98nBMSQk2ZQr84Afw+OOGfktEXZ2zDbAc\nGACsA+YAQ4Bl9dr0B5YCHxL+kagC+jWyL1fnlLRHEybA1VeHryefHHdv4tXS1Tmjjvj7AiuB1cAO\nYBwwuEGb1wmhDzALODziMSUl1KOPwnXXwQsvGPpRRK3xHwasqfd+LXDCbtpfATwX8ZiSEuiuu8LT\ns6ZPh27d4u5NaYsa/M2pzZwGXA6c1FSDqqqqz15XVlZS6XJ6UuKl0/Czn4W19F95BTp1irtH8Uml\nUqRSqcj7iVrj70eo2Q/Mvh8BpIHbGrTrBUzItlvZxL6s8Uv6nO3b4eKL4f33YeJEOPjguHtUXOKq\n8c8FjgY6AXsBFwKTG7T5Z0Lo/5imQ1+SPucf/wjr7rRrB1OnGvr5FDX4dwLDgSmEmTtPEWb0DMtu\nAL8EDgLuBRYAsyMeU1KZq6mBk06CU06BJ5+EvfeOu0flxYetSyoqs2fD4MHwi1/AtdfG3Zvi5sPW\nJZW8CRNg2DB46CH4/vfj7k35MvglxS6dDsspP/RQmL3z7W/H3aPyZvBLitWWLXDJJfD3v4cyz6GH\nxt2j8ueyzJJiU1MD/ftDhw7hxixDvzAMfkmxePFFOPFEuOYaeOABZ+4UkqUeSQWVycDdd8Mtt8BT\nT4WHqKiwDH5JBbNlC1x1VXhq1syZyV5+IU6WeiQVxKJFcPzxsP/+8Prrhn6cDH5JrSqTgTFjwvIL\nI0eGev4++8Tdq2Sz1COp1Xz0UXhoyqJFYWXNrl3j7pHAEb+kVrJ4cbgRa++9YdYsQ7+YGPyS8iqd\nhtGj4fTT4cYbYexY2HffuHul+iz1SMqbdevg8sth0yZ49VU45pi4e6TGOOKXlBfjx0Pv3mE5ZUO/\nuDnilxTJpk0wfDjMmQPPPgt9+sTdI+2JI35JLTZ9OvTqBe3bw4IFhn6pcMQvqdk2bYKf/zwsoTxm\nDAwcuOefUfFwxC+pWSZMgB49oG1bWLLE0C9Fjvgl5WTdulDLr64Oi6udfHLcPVJLOeKXtFvpNNx/\nPxx3HPTsGWr5hn5pc8QvqUkLF4ZR/s6d4UJuz55x90j54Ihf0hds2BAekHL22eGxiK++auiXE4Nf\n0md27YLf/x66dQsXb6urw/r5bdrE3TPlk6UeSUBYPfP668Oc/BdfDPPzVZ4MfinhVq6Em24KD0e5\n4w744Q+hoiLuXqk1WeqREmr9erjuOujXD775zfA4xB/9yNBPAoNfSpjNm8OTsLp3D2vlV1eH5ZP3\n2y/unqlQDH4pIT75BP7rv8KqmatWwbx58NvfQocOcfdMhZaP4B8IVAMrgBsa+bwr8DrwMfB/8nA8\nSc3w8cdwzz3QpUtYW2fKFHjsMR92nmRRL+62AUYDA4B1wBxgMrCsXpsNwPXAuRGPJakZtm8PDza/\n/fawTv6f/gR9+8bdKxWDqCP+vsBKYDWwAxgHDG7Q5n1gbvZzSa1s61a480446ihIpWDyZHjmGUNf\ndaKO+A8D1tR7vxY4IeI+JbXA++/DvfeGG7BOOQVeeCHM1pEaihr8mbz0Iquqquqz15WVlVRWVuZz\n91JZWr4cfve7sGLm+eeHNXW6d4+7V2oNqVSKVCoVeT9RZ+z2A6oIF3gBRgBp4LZG2o4EPgLubGJf\nmUwmr/+OSGUrk4G//CXMypk5M6yrc+21cMghcfdMhVQRbrpodo5HHfHPBY4GOgHvABcCQ5po620h\nUkTbtsG4cWGWzkcfwb/9Wxjp77NP3D1TKclHGA8CRhFm+IwFbgGGZT+7HziUMNvnAMJvA1uA7oTR\nf32O+KUmLF0a1sR/4gk48US4+moYNAi+5J04idbSEX8xjcINfqmeTz4Jjzm87z546y248sqwHXFE\n3D1TsTD4pTKQycD8+fDoo6Gk06tXGN0PHgzt2sXdOxWbuGr8kvJg3Tp48skQ+Nu3h4efvP46dO4c\nd89UjhzxSzHZvDncXPX44zB7dpiKeemlcNJJ1u6VG0s9UgnYvDncRTt+PMyYAaeeCkOHhlLOvvvG\n3TuVGoNfKlK1Yf/HP4abq049NTzs5JxzwtOupJYy+KUismpVCPvJk2HWrBD2P/qRYa/8MvilGO3a\nFQL+mWfC9v778N3vwve/D2eeCV/5Stw9VDky+KUCW7UKpk6FadNCCeeww0LQn3NOWAnTC7RqbQa/\n1Mo2bgzr40ybFrbNm8No/swzYcCAEPxSIRn8Up6tXw8vv1y3rVoF/fvXhX3Pno7qFS+DX4ognQ7L\nIsycCa+9FoL+vffg5JPDhdnvfCc8xcq7Z1VMDH6pGTZtChdjZ84M26xZcOCB0K9f2L7znbBcQps2\ncfdUaprBLzVh48aw/s28eWGbPz+M5r/97bqg79cPDj007p5KzWPwK/HS6VCHf+ONum3+fNiwAY47\nDo4/Pmy9e8OxxzqaV+kz+JUYmUwYsS9dGrYlS0LIL1kCBx8cSjQ9e4atd284+mgvwqo8GfwqOzt3\nwurV4aJrdXVd0C9bFi6ydu8O3bpBjx7hoeI9e3pXrJLF4FdJSqfhnXdg5cqwLV8egv6tt0LZ5mtf\ng2OOCaWZ7t3rwv6f/inunkvxM/hVtD75BP72N/if/wlbTU1d0K9aFWbTdO4MXbqEgK8N+s6dfZas\ntDsGv2KTTsO774YQX706fF21qi7o16+Hww+Ho46CI48MAd+lSwj2zp1dx0ZqKYNfrSadDhdTV68O\n29/+Vhfwq1fD22/DQQdBp04h2Dt1CiFfux1+OLT1WW9S3hn8arF0OozKa4O9NtBrt7ffDhdNO3UK\nD/qu/XrkkWE74ghLMlIcDH41KZMJNzHVlmAabm+/DQccEAK9dqsduR9xRNh8OpRUfAz+hNu1C9as\nCRdOa7faC6k1NaFNbY29NtTrv95vvzh7L6klDP4ESKdh7dow1XHFijArZsWKsK1eDR06hIumRx0V\nLprWfu3cOdTgK4rpT1tSZAZ/Gdm+PYT78uXhxqXabfnyMPXxmGPC3aj1N6c+Sslj8JegTz8NAb9k\nCbz5Zt3XNWtCkHftWrcde2zYDjgg7l5LKhYGf5HbsAEWLYIFC2DhwrCtXBkunH7jG2HZgW98I2xd\nurjuu6Q9M/iLyN//DnPmhG3+/BDyH34Y1pM57ri6rXt3+PKX4+6tpFIVZ/APBEYBbYAHgdsaaXM3\nMAjYBlwGLGikTUkG/0cf1YV87bZpU1jrvU+f8PVb3wozZ1whUlI+xRX8bYDlwABgHTAHGAIsq9fm\nX4Dh2a8nAHcB/RrZV0kE/9q18OqrdVt1dRjJ9+0bgr5Pn1CqMeQltbaWBn/UG+n7AiuB1dn344DB\nfD74zwEezb6eBbQHDgHWRzx2q8tkwlz4GTPC9te/wtatcNJJYbv77vBgD8s1kkpJ1OA/DFhT7/1a\nwqh+T20Op0iDf82aEPLTp4evO3bA6afDaafBL38ZplI6H15SKYsa/LnWZhpGZdHUdLZvh1QKnn8e\nXngBPvgghPxpp8GIEQa9pPITtRK9DuhY731Hwoh+d20Oz37vCyoqquptKSoqoKqq8QNXVYVAbrg1\nt3379nDLLeGBH089FRYrGz8+fO3aNdTqo+zf9ra3ve3z1f6yy1JUVVV9trVURYt/MmhLuLh7BvAO\nMJvdX9ztR5gBVNCLu+k0zJoFEyfCpEmhTj9wIAwaBAMG+Lg+SaUprou7OwmhPoUww2csIfSHZT+/\nH3iOEPorga3ATyIeMyeffhpq9BMnwp//DF/9KvzgB/CHP4QHcFdE/SdPkkpUMcVf5BH/jh3w0kvw\n3/8Nzz4bSjXnnhsC/+ij89RLSSoSib1zN52G114LI/k//jGscTNkCFxwAXz9663QS0kqEnGVemLz\n5pvw2GMwbhzsvz8MHQozZ4aliCVJTSup4N+8Ocy8GTs2zLe/+OJQ0unZM+6eSVLpKPpSTyYTlkYY\nOzbMyDntNLjiCjj7bB/gLSnZyq7G/+GH8PDDcN99YQbOFVeEEf4hh8TYQ0kqImVT41+6FEaPDrX7\ns8+GBx8M6+I4/VKS8qOogn/AgHDR9qqrwtOonJUjSflXTOPozBNPZLjgAth777i7IknFr+xq/JKk\n3Wtp8Pu4EElKGINfkhLG4JekhDH4JSlhDH5JShiDX5ISxuCXpIQx+CUpYQx+SUoYg1+SEsbgl6SE\nMfglKWEMfklKGINfkhLG4JekhDH4JSlhDH5JShiDX5ISxuCXpISJEvwHA9OAt4CpQPsm2j0ErAcW\nRziWJClPogT/vxOC/xjgpez7xjwMDIxwnERJpVJxd6FoeC7qeC7qeC6iixL85wCPZl8/CpzbRLtX\ngA8iHCdR/Etdx3NRx3NRx3MRXZTgP4RQwiH79ZDo3ZEktba2e/h8GnBoI9+/qcH7THaTJBW5igg/\nWw1UAu8BXwNmAF2baNsJeAbouZv9rQQ6R+iPJCVNDdCluT+0pxH/7kwGLgVuy36dFGFf0ILOS5IK\n62DgRb44nfPrwP+r1+4PwDvAJ8Aa4CcF7KMkSZKkQhpIuDawArihiTZ3Zz9fBHyrQP2Kw57OxUWE\nc/AG8CrQq3BdK7hc/l4A9AF2AucVolMxyeVcVAILgCVAqiC9iseezkUH4AVgIeFcXFawnhVWLjfB\nFm1utiFcwO0EtCP8YXVr0OZfgOeyr08AZhaqcwWWy7noDxyYfT2QZJ+L2nbTgWeB8wvVuQLL5Vy0\nB94EDs++71CozhVYLueiCrgl+7oDsIFo1y2L1SmEMG8q+Judm4Vcq6cv4Q9yNbADGAcMbtCm/k1h\nswh/ycvx/oBczsXrwIfZ17Oo+x+93ORyLgCuB54G3i9Yzwovl3MxFPgTsDb7/h+F6lyB5XIu3gUO\nyL4+gBD8OwvUv0La002wzc7NQgb/YYSLu7XWZr+3pzblGHi5nIv6rqDuX/Ryk+vfi8HAvdn35XrP\nSC7n4mjCxIoZwFzg4sJ0reByORdjgB6EySOLgP9dmK4VnWbnZiF/Lcr1f9aG9xaU4//kzflvOg24\nHDiplfoSt1zOxSjCWlAZwt+PKPefFLNczkU7oDdwBrAv4TfDmYT6bjnJ5VzcSCgBVRLuAZoGfBPY\n0nrdKlrNys1CBv86oGO99x2p+3W1qTaHZ79XbnI5FxAu6I4h1PjLdb2jXM7F8YRf9SHUcgcRfv2f\n3Oq9K6xczsUaQnlne3Z7mRB25Rb8uZyLE4H/m31dA6wCjiX8JpQkRZ2bbQl/OJ2Avdjzxd1+lO8F\nzVzOxT8Tapz9CtqzwsvlXNT3MOU7qyeXc9GVcP9MG8KIfzHQvXBdLJhczsVvgZHZ14cQ/mE4uED9\nK7RO5HZxtyhzcxCwnBBoI7LfG5bdao3Ofr6I8CttudrTuXiQcLFqQXabXegOFlAufy9qlXPwQ27n\n4meEmT2Lgf9V0N4V1p7ORQfCUjCLCOdiaKE7WCC1N8F+SviN73KSm5uSJEmSJEmSJEmSJEmSJEmS\nJEmSVNz+P4qdeORFAjw6AAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x85a6910>"
       ]
      }
     ],
     "prompt_number": 70
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Z = np.roots([1,B-1,A-2*B-3*B*B,-A*B+B**2+B**3])\n",
      "rho = p/(Z*R*T) # Z = p/(rho*R*T)\n",
      "rho"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 71,
       "text": [
        "array([  69.42536782-63.64922759j,   69.42536782+63.64922759j,\n",
        "        518.88328530 +0.j        ])"
       ]
      }
     ],
     "prompt_number": 71
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "rho_EOS"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 72,
       "text": [
        "495.5913139092507"
       ]
      }
     ],
     "prompt_number": 72
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Mixtures"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "components = ['Ethanol','Water']\n",
      "z = [0.4, 0.6]\n",
      "#T = 562.8969961 #[K]\n",
      "T = 400 #[K]\n",
      "p = 12e6 #[Pa]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 117
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "omega = np.array([Props(f,'accentric') for f in components]) #[-]\n",
      "R = 8.314472 #[J/mol/K]\n",
      "Tc = [Props(f,'Tcrit') for f in components] #[K]\n",
      "pc = [Props(f,'pcrit')*1000 for f in components] #[Pa]\n",
      "\n",
      "A = 1.1\n",
      "B = 0.25\n",
      "C = 0.2\n",
      "D = 1.2\n",
      "E = 0.01\n",
      "\n",
      "kappa_0 = 0.378893 + 1.4897153*omega-0.17131848*omega**2+0.0196554*omega**3\n",
      "kappa_1 = [kappa_dict[f] for f in components]\n",
      "\n",
      "a,b = 0,0\n",
      "for i in range(len(components)):\n",
      "    \n",
      "    b += z[i]*0.077796*R*Tc[i]/pc[i]\n",
      "    \n",
      "    for j in range(len(components)):\n",
      "        kappa_i = kappa_0[i]+kappa_1[i]*(sqrt((A-D*(T/Tc[i]+B))**2+E)+A-D*(T/Tc[i]+B) )*sqrt(T/Tc[i]+C)\n",
      "        alpha_i = (1+kappa_i*(1-sqrt(T/Tc[i])))**2\n",
      "        a_i = 0.457235*R**2*Tc[i]**2/pc[i]*alpha_i\n",
      "        \n",
      "        kappa_j = kappa_0[j]+kappa_1[j]*(sqrt((A-D*(T/Tc[j]+B))**2+E)+A-D*(T/Tc[j]+B) )*sqrt(T/Tc[j]+C)\n",
      "        alpha_j = (1+kappa_j*(1-sqrt(T/Tc[j])))**2\n",
      "        a_j = 0.457235*R**2*Tc[j]**2/pc[j]*alpha_j\n",
      "        \n",
      "        if i != j:\n",
      "            k_ij = -0.3\n",
      "        else:\n",
      "            k_ij = 0\n",
      "            \n",
      "        a += z[i]*z[j]*(a_i*a_j)**0.5*(1-k_ij)\n",
      "\n",
      "AA = a*p/(R**2*T**2)\n",
      "BB = b*p/(R*T)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 118
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "coeffs = [1,BB-1,AA-2*BB-3*BB**2,-AA*BB+BB**2+BB**3]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 119
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Z = np.linspace(0,1)\n",
      "def f(Z):\n",
      "    return sum([a*Z**b for (a,b) in zip(coeffs,reversed(range(len(coeffs))))])\n",
      "plt.plot(Z,f(Z))\n",
      "plt.axhline(0,linestyle = '--')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 120,
       "text": [
        "<matplotlib.lines.Line2D at 0x97cd2d0>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEACAYAAAC08h1NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHK5JREFUeJzt3Xu8jXXax/HPzpBUkjQKGYVE0oEMnaxeowaFRk0lHZ5S\nVMxUzzOFmhlrpoNDj0PKIKeohnSYyamEWh4ktRuknMXkULscJsRmH9bzx7Wwbdvey77Xvk/r+369\n9qu19rr3ui93XPu3rt/vvn4gIiIiIiIiIiIiIiIiIiIiIiIiIiIBNw7IApYf4/UuwDLgC2Ah0MSl\nuEREpIxcDVzKsRN/S+C0xOM2wCduBCUiImWrDsdO/AWdDmwu21BERKQ4J7h8vq7ATJfPKSIiZaAO\nJY/4rwVWYKN+ERHxyM9cOk8TYDRW499Z+MW6devG169f71IoIiKhsR6od7w/5EappzbwDnAnsK6o\nA9avX088HtdXPE7fvn09j8EvX7oWuha6FsV/AXVLk5RTMeKfBLQCqgGbgL5A+cRro4A/Y+WdEYnv\n5QDNU3BeEREphVQk/s4lvH5/4ktERHzA7VU9UoJIJOJ1CL6ha3GYrsVhuhbOZXgdQEI8Ua8SEZEk\nZWRkQCnyuEb8IiJpRolfRCSAcnJK/7NK/CIiAROPw0MPlf7nlfhFRALm6adhyZLS/7wSv4hIgIwf\nD6+8AjNmlP49tKpHRCQgZs2Ce+6BefOgQYPSr+pxq1ePiIg48K9/wV13wT//aUnfCZV6RER8buNG\naN8eRo6EK65w/n5K/CIiPrZjB7RtC717Q6dOqXlP1fhFRHxq3z64/npo0QKef/7o10tb41fiFxHx\nodxcuOUWqFQJXnsNTiiiPqPJXRGRkIjHoXt3yM6GKVOKTvpOKPGLiPhMnz7w1Vcwdy5UqJD691fi\nFxHxkUGDYOpUmD8fTj65bM6hxC8i4hMTJsCwYbBgAZxxRtmdR5O7IiI+MG0adOsGH30EF1yQ3M9o\ncldEJKDmz4euXa3/TrJJ3wndwCUi4qElS2zZ5uuvw+WXu3NOJX4REY+sXAnt2sGIEXDdde6dV4lf\nRMQDX39td+UOHJi6VgzJUuIXEXHZ5s3QujU8+aR13HSb08Q/DsgClhdzzDBgLbAMuNTh+UREAu37\n7y3pP/SQs+0TnXCa+McDbYp5vR1QD6gPdANGODyfiEhg7dxp5Z1bb4XHH/cuDqeJfz6ws5jXOwAT\nEo8XA1WA6g7PKSISOLt320TutdfCX/7ibSxlXeOvCWwq8HwzUKuMzyki4it790LHjtC4MQweDBke\n3zrrxg1chf+IRd6iG41GDz2ORCJEIpGyi0hExCX79lnSr1XLdtBykvRjsRixWMxxTKn4vVMHmAZc\nVMRrI4EYMDnxfBXQCpsQLkgtG0QkdLKz4aabrO/OxIlQrlxq37+0LRvKutQzFbg78bgF8B+OTvoi\nIqGzfz/cfDOcdpo1X0t10nfCaalnEjaCr4bV8vsC5ROvjQJmYit71gE/Afc6PJ+IiO8dOAC//S2c\ndJLtnvUzn3VFU3dOEZEUysmx5Zpgu2eVL1/88U6oO6eIiMdycqBzZ8jLg7feKtuk74QSv4hICuTk\nQJcutnTzH/8omy0TU0WJX0TEoQMH4Pbb7b/vvAMnnuh1RMVTkzYREQcOrt7Jz7ekX7Gi1xGVTIlf\nRKSU9u2zdfoVK8Kbb/q7vFOQEr+ISCns3QsdOkCVKjBpkn8ncouixC8icpz27IEbboCzz4ZXX/Xf\nOv2SKPGLiByH3buhbVs491wYPz54SR+U+EVEkrZjh22i0qgRjBnjrzYMx0OJX0QkCd99B5EIXH21\nddk8IcDZM8Chi4i445tv4JprrP/O889730/fKSV+EZFirFljo/wePeBPfwp+0gfduSsickzLltlE\n7jPPwH33eR1N6ijxi4gU4ZNPbOesF1883G0zLJT4RUQKmTvXumy+8optkB42qvGLiBQwZQrccYe1\nVQ5j0geN+EVEDhk+HPr1g9mzoUkTr6MpO0r8IpL24nHo2xcmT4b58+2u3DBT4heRtJaXBw8/DJ9/\nDgsWwM9/7nVEZU+JX0TSVna27Zq1axd89BGceqrXEblDk7sikpb+8x9bo1++PEyfnj5JH5T4RSQN\nbdoEV11lE7h//7v/t0pMNSV+EUkrS5fCFVdA167wwgvBbrZWWqn4I7cBVgFrgV5FvF4NeB9YCnwJ\n/FcKzikictxmzYLrr4chQ+Cxx7yOxjtO2w2VA1YDrYEtwGdAZ2BlgWOiwIlAH+yXwGqgOpBb4Jh4\nPB53GIqIyLGNGwdPPglvvw1XXul1NKmRYR3jjjuPO13V0xxYB2xMPJ8MdOTIxP8tcPBWiMrAdo5M\n+iIiZSYeh2gUXnsN5s2DBg28jsh7ThN/TWBTgeebgV8WOmY08CGwFTgVCFm7IxHxqwMHoFs3WLEC\nPv4Yqlf3OiJ/cJr4k6nPPInV9yNAXWA2cDGwu+BB0Wj00ONIJEIkEnEYmoiks+3boVMnOP10W6N/\n8sleR+RcLBYjFos5fh+nNf4WWA2/TeJ5HyAfGFDgmJnAs8DCxPO52CRwZoFjVOMXkZRZvRpuvNES\nf79+4V25U9oav9PLkQnUB+oAFYDbgKmFjlmFTf6CTeo2AL52eF4RkSJ9+KFtk9i7NwwYEN6k74TT\nUk8u0BOYha3wGYtN7HZPvD4KeA4YDyzDftE8AexweF4RkaOMGQNPPWXN1q691uto/Msvu0eq1CMi\npZaXZyP8d9+19gvnn+91RO7wajmniIindu2Cu+6CH3+07RKrVvU6Iv9T9UtEAmv9emjZEs46Cz74\nQEk/WUr8IhJIc+ZYz52ePWHUKKhQweuIgkOlHhEJlHjcmqsNGGD747Zq5XVEwaPELyKBkZ0NDz5o\nHTYXLYI6dbyOKJhU6hGRQPj2W4hEYO9eWLhQSd8JJX4R8b2FC6FZM7sb9403wtF+wUsq9YiIb8Xj\n8NJL8MwzMH48tGvndUThoMQvIr60dy907w7Ll1tnzbp1vY4oPFTqERHfObg+H5T0y4ISv4j4ysyZ\ntj7/gQdg4kSoVMnriMJHpR4R8YW8PPjrX2HsWHjnnfBsj+hHSvwi4rmsLOjSBfLzITPTWjBI2VGp\nR0Q8NW8eNG1q5Z3Zs5X03aARv4h4Ij/f2i688AJMmAC//rXXEaUPJX4Rcd327XD33dZKOTMTatXy\nOqL0olKPiLjq44/hssugUSPbBF1J330a8YuIK/LyoH9/GDYMRo+GDh28jih9KfGLSJnbuhXuvNPq\n+p9/rlG+11TqEZEyNWOGlXYiEZg7V0nfDzTiF5EysX+/bYD+9tvw5ptw9dVeRyQHKfGLSMqtWmU3\nZNWubZumaC9cf1GpR0RSJh6HkSPhqqvg/vut9YKSvv+kYsTfBhgKlAPGAAOKOCYCDAHKA9sSz0Uk\nRL7/Hrp2tYncBQvgggu8jkiOxemIvxzwEpb8GwGdgYaFjqkCDAfaA42BWxyeU0R8Zvp0uPhiuOgi\n2wtXSd/fnI74mwPrgI2J55OBjsDKAsfcAbwNbE483+bwnCLiE3v3wh/+YK2Up0zRBG5QOB3x1wQ2\nFXi+OfG9guoDVYGPgEzgLofnFBEf+PRTW6a5ezcsW6akHyROR/zxJI4pD1wG/AqoBCwCPgHWFjwo\nGo0eehyJRIhEIg5DE5GycOCA9c0fPRpefBFuvdXriNJHLBYjFos5fp8Mhz/fAohiNX6APkA+R07w\n9gJOShwHNgH8PvBWgWPi8Xgyv0NExEvLlsE998A551jiVwtlb2VkZEAp8rjTUk8mVsqpA1QAbgOm\nFjrmXeAqbCK4EvBLYIXD84qIi3Jz4bnnoHVrePRRmDpVST/InJZ6coGewCwssY/FJna7J14fBazC\nRvhfYJ8GRqPELxIYq1fbKP/kk63PTu3aXkckTjkt9aSKSj0iPpObC0OG2GYp0Sg8/DCcoFs+faW0\npR61bBCRo3z1Fdx7L5xyiq3eOe88ryOSVNLvbxE5JCcHnnnGOml27Qpz5ijph5FG/CICWDO1e++F\n6tVVyw87jfhF0lx2Nvzxj3D99bZi5733lPTDTiN+kTQ2bx5062Y9dpYuhRo1vI5I3KDEL5KGdu6E\nxx+HWbPgpZegY0evIxI3qdQjkkbicXjjDbjwQqhY0VbvKOmnH434RdLEv/8NPXrAxo22HWLLll5H\nJF7RiF8k5HJyYOBAaNoUWrSAf/1LST/dacQvEmLz58NDD1lTtcWLoW5dryMSP1DiFwmhH36AJ56w\nG7CGDIGbb4YMvzRoEc+p1CMSIvn5MGYMNG4Mp58OK1bALbco6cuRNOIXCYnMTOjZ05L8rFlwySVe\nRyR+pRG/SMBt2wbdu8ONN9p/Fy5U0pfiKfGLBFReHowcCY0a2Zr8Vaus145aJ0tJVOoRCaBFi2xN\n/imn2ARukyZeRyRBosQvEiBbtkDv3vDhh/D889C5syZu5fjpQ6FIAOzbB08/bSP72rVtO8Q77lDS\nl9LRiF/Ex+JxePNNW5PfrJmt3Dn3XK+jkqBT4hfxqSVL4JFHYNcueOUV2xVLJBVU6hHxmS1bbHVO\n27Zw5522G5aSvqSSEr+IT+zZA337Wh3/rLNgzRrbJKVcOa8jk7BR4hfxWF4ejB0LDRrAunXWPbNf\nP6hc2evIJKxSkfjbAKuAtUCvYo67HMgFOqXgnCKhMHs2XHaZ1fD/8Q94/XX4xS+8jkrCzunkbjng\nJaA1sAX4DJgKrCziuAHA+4AWoEnaW7IEevWyTVH694ff/EZLM8U9Tkf8zYF1wEYgB5gMFLWR2++A\nt4AfHJ5PJNA2brQJ23btLNl/9RV06qSkL+5ymvhrApsKPN+c+F7hYzoCIxLP4w7PKRI427bBY4/Z\nWvz69WHtWtsgpXx5ryOTdOS01JNMEh8K9E4cm8ExSj3RaPTQ40gkQkTr1yQEfvoJXngBBg+G226z\nEX716l5HJUEVi8WIxWKO38fpB8wWQBSb4AXoA+Rj9fyDvi5wnmrAXuABbC7goHg8rg8CEh4HDsDL\nL8Nzz8E111i7hfr1vY5KwibDaoTHncedjvgzgfpAHWArcBvQudAx5xV4PB6YxpFJXyQ08vJsZU7f\nvtCwIcyYAZde6nVUIkdymvhzgZ7ALGzlzlhsRU/3xOujHL6/SCDE4/Duu/DHP0KVKjBxIlx9tddR\niRTNL2sJVOqRQIrH4YMP4M9/huxsK+20a6dVOuIOr0o9ImkrFoM//clW7ESj8NvfavcrCQYlfpHj\ntGiRJfwNGyzh33GH+ulIsGh8IpKkzEy44QZblnn77bbH7V13KelL8Cjxi5QgMxPat4ebbrJWyWvX\nwv336+YrCS4lfpFjKJjwf/1r65zZsyeceKLXkYk4o8QvUsjnn0OHDtCxI1x//eGEX7Gi15GJpIYS\nv0jC4sVw442W9K+7Dtavh9/9TglfwkereiTtLVhgLRVWrrRWyW+9pWQv4abEL2kpHrd1+H/9K/z7\n39CnD9xzD1So4HVkImVPiV/SSjwO778Pzz4LWVnw1FPQpYtW6Eh6UeKXtJCfb1sbPvcc7N8PTz4J\nt94KP9O/AElD+msvoZaTA5Mm2eblp55qPXXat1drBUlvSvwSStnZtoH5gAG2efmwYdC6tZqniYAS\nv4TMjz/CyJEwdChcdhm89hpceaXXUYn4ixK/hEJWlm1x+PLL0KYNzJoFTZp4HZWIP6nSKYG2YQP0\n6GG7Xf34I3z6qY3ylfRFjk2JXwJpyRLo3BmaNYPKle3mq+HD4bzzSv5ZkXSnxC+BEY/DnDnWP6d9\ne2ja1Eb8/fpB9epeRycSHKrxi+/l5sLbb8PAgbBvHzzxhG1+ortsRUrHL4vbtOeuHOWnn2DcOBgy\nBGrUsD46N9ygNfgiB2nPXQmNrCx48UUYNQquuQZefx1atvQ6KpHw0NhJfGP1aujWDS64ALZvh48/\nthKPkr5IamnEL56Kx2H+fBg0yDYxf/hhWLMGzjzT68hEwisVI/42wCpgLdCriNe7AMuAL4CFgFZY\nC7m58MYb0Ly57V/bti1s3AjRqJK+SFlzOrlbDlgNtAa2AJ8BnYGVBY5pCawAfsR+SUSBFoXeR5O7\naWL3bhgzxu6yrV0b/ud/1DRNpLS8mtxtDqwDNiaeTwY6cmTiX1Tg8WKglsNzSgB9841N2I4bZ83S\npkyx0b6IuM/pOKsmsKnA882J7x1LV2Cmw3NKgGRm2pr7Sy+FvDzbyPxgiUdEvOF0xH889ZlrgfuA\nInslRqPRQ48jkQiRSMRJXOKhvDyYNg0GD7ZtDR95BEaMgNNO8zoykWCLxWLEYjHH7+O0xt8Cq9m3\nSTzvA+QDAwod1wR4J3HcuiLeRzX+ENizx3rgDx0KZ5xh9ftOnbTLlUhZ8arGnwnUB+oAW4HbsMnd\ngmpjSf9Oik76EnCbNh2u30ciMGECXHGFNj0R8SuniT8X6AnMwlb4jMUmdrsnXh8F/Bk4HRiR+F4O\nNiksAffpp9ZOYdYsuOce+OwzOPdcr6MSkZL4ZUymUk9A5OXBP/9pCX/zZvj976FrV9XvRbygXj1S\npnbtgrFjbe/aGjXg0UfhN79R/V4kiPTPVoq1YYMl+4kTrQ++lmKKBJ/ul5SjxOOwYAHcfDNcfrn1\nvV+6FCZNUtIXCQON+OWQAwfgzTetfr9rl62/nzABTjnF68hEJJU0uSts326974cPt5bIjz6qDU9E\ngqC0k7v6p53GVq2CBx+EevVg7VqYORPmzlXTNJGwU6knzRzcsHzIEOub8+CDsHIlnHWW15GJiFuU\n+NNEdrZtYTh0qD1/9FF45x2oWNHbuETEfUr8IZeVBX/7G4wcCc2a2Uj/V79SOwWRdKZKbkh98QXc\ne69N1n7/PcybBzNmWC98JX2R9KYRf4jk58N779mofuVK6NED1q2zTpkiIgcp8YfA3r3w6quW8CtV\ngv/+b7j1VrvxSkSkMCX+APvuO1t7P2oUtGxpdfxWrVTKEZHiqcYfQAfr940awY4d1l7h3XetF76S\nvoiURCP+gIjHre/9oEGwYgX07Gn1+6pVvY5MRIJGid/n9u+39feDB0O5crad4e23q34vIqWnxO9T\n27fbBuXDh8Mll9iNV1p/LyKpoBq/z6xbZ8sw69WzXvhz5tgSTa2/F5FUUeL3iU8+sf73LVtClSq2\nDn/sWLjwQq8jE5GwUanHQ/n5MG0aPP88bNli6+/V/15EypoSvweys20rw0GDoHJlePxx6NRJ+9eK\niDuUaly0c6dN2L74IjRtajde6YYrEXGbavwu+OYbK+PUrQtr1sDs2TB9um64EhFvpCLxtwFWAWuB\nXsc4Zlji9WXApSk4ZyB8+SXcfbctxzzhBFi2DF55BRo39joyEUlnThN/OeAlLPk3AjoDDQsd0w6o\nB9QHugEjHJ7T9xYutO0Lr7sOGjaEr7+G//1fOOccryMTEXFe428OrAM2Jp5PBjoCKwsc0wGYkHi8\nGKgCVAeyHJ7bV+JxW2/frx9s3WoTtlOmwEkneR2ZiMiRnCb+msCmAs83A79M4phahCTx5+Zagu/f\n38o5vXvDLbdohY6I+JfTUk88yeMKT2Ee9XMZGdECXzEyMiAaLfrNolGbFC385ebx+/fDyy9DgwbQ\npw8sX241/M6doXx59+PR8Tpex4f/+FgsRjQaPfRVWhml/knTAohiNX6APkA+MKDAMSOBGFYGApsI\nbsWRI/54PJ7s7xBv7dljCX/QIJu07dMHrrrK66hEJB1lZGRAKfK40xF/JjZpWweoANwGTC10zFTg\n7sTjFsB/CGCZZ+dOePppOO88a68wfbrtYaukLyJB47QSnQv0BGZhK3zGYhO73ROvjwJmYit71gE/\nAfc6PKerfvjBWiK//DJ07Ajz51t5R0QkqJyWelLFd6Web7+1JZjjx1v/+1694Be/8DoqEZHDvCr1\nhM6mTba71YUXQl6eTdr+7W9K+iISHkr8CRs2QLducPHFUKmStUUeOhRq1vQ6MhGR1Er7xL9hAzzw\nADRrBmeeab10Bg6E6tW9jkxEpGykbeLfsAHuv98SfvXqlvCffRaqVfM6MhGRspV2ib9gwj/7bFi7\nFp55Bs44w+vIRETckTaJ/5tvrIZfMOE//TRUrep1ZCIi7gp94t+61VbpXHKJjerXrFHCF5H0FtrE\nn5Vlm580bgwnngirVlnnTJV0RCTdhS7x79hh/XMaNoScHNsMZdAg+PnPvY5MRMQfQpP49+yxVTnn\nnw/btsHSpba3bY0aXkcmIuIvgU/8+/fDsGFQr57dZfvxxzB6NNSu7XVkIiL+FNjtQnJz4dVXrUf1\nRRfB++/bBK6IiBQvcIk/HoepU62Of+aZ8Pe/w5VXeh2ViEhwBCrxL1wITzwBu3db58y2bW1nGhER\nSV4gEv+KFTbCX7rU1uB36QLlynkdlYhIMPl6cnfLFmuvEIlAq1awejXcfbeSvoiIE75M/Hv2QN++\n0KSJNU1bs8ZuxqpY0evIRESCz1eJPzfXlmKefz58/TUsWQL9+0OVKl5HJiISHr6p8b/3Hjz+uI3w\np061ZmoiIpJ6flkTE2/QIM7AgdC+vVbqiIgko7R77volxcYPHIhTvrzXYYiIBEfgE388Hvc6BhGR\nQClt4ncyuVsVmA2sAT4AipqCPQf4CPgK+BL4vYPziYhICjhJ/L2xxH8+MDfxvLAc4DHgQqAF0ANo\n6OCcoReLxbwOwTd0LQ7TtThM18I5J4m/AzAh8XgCcFMRx3wHLE083gOsBNQouRj6S32YrsVhuhaH\n6Vo45yTxVweyEo+zEs+LUwe4FFjs4JwiIuJQSev4ZwNnFfH9pwo9jye+juUU4C3gEWzkLyIiHnGy\nqmcVEMHKOWdjk7gXFHFceWA68B4w9BjvtQ6o6yAWEZF0tB6o5+YJBwK9Eo97A/2LOCYDmAgMcSso\nEREpO1WBORy9nLMGMCPx+CogH5vgXZL4auNumCIiIiIi4po22NzAWg6XiQoblnh9GbYKKKxKuhZd\nsGvwBbAQaOJeaK5L5u8FwOVALtDJjaA8ksy1iGCfnr8EYq5E5Y2SrkU14H2sovAl8F+uReaucdjK\nyeXFHOPbvFkOm8Stg034LuXom7naATMTj38JfOJWcC5L5lq0BE5LPG5Del+Lg8d9iC0UuNmt4FyW\nzLWogt0JXyvxvJpbwbksmWsRBfolHlcDtuOjjsMpdDWWzI+V+I87b7rZj7859j9yI3ZH72SgY6Fj\nCt4Uthj7S17S/QFBlMy1WAT8mHi8mMP/0MMmmWsB8DtsSfAPrkXmvmSuxR3A28DmxPNtbgXnsmSu\nxbdA5cTjyljiz3UpPjfNB3YW8/px5003E39NYFOB55sT3yvpmDAmvGSuRUFdOfwbPWyS/XvRERiR\neB7Wjn7JXIv62MKKj4BM4C53QnNdMtdiNNYOZitW4njEndB857jzppsfi5L9x1r43oIw/iM/nj/T\ntcB9wJVlFIvXkrkWQ7Elw3Hs74dfusqmWjLXojxwGfAroBL2yfATrL4bJslciyexElAEuw9oNnAx\nsLvswvKt48qbbib+LVi3zoPO4fDH1WMdUyvxvbBJ5lqATeiOxmr8xX3UC7JkrkVT7KM+WC23Lfbx\nf2qZR+euZK7FJqy8sy/x9X9Ysgtb4k/mWlwBPJt4vB7YADTAPgmlE1/nzZ9h/3PqABUoeXK3BeGd\n0EzmWtTGapwtXI3Mfclci4LGE95VPclciwuw+2fKYSP+5UAj90J0TTLXYjDQN/G4OvaLoapL8bmt\nDslN7voyb7YFVmMJrU/ie90TXwe9lHh9GfaRNqxKuhZjsMmqgze+fep2gC5K5u/FQWFO/JDctfgD\ntrJnOeHe46Kka1ENmIbliuXYxHcYTcLmMQ5gn/juI33zpoiIiIiIiIiIiIiIiIiIiIiIiIiIiIiI\niIiIv/0/GIpr6V1/g5sAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x966b4d0>"
       ]
      }
     ],
     "prompt_number": 120
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Z = np.roots(coeffs)\n",
      "rho = p/(Z*R*T) # Z = p/(rho*R*T), rho in [mol/m^3]\n",
      "rho"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 121,
       "text": [
        "array([  1234.16865590-3249.94868206j,   1234.16865590+3249.94868206j,\n",
        "        24827.37270034   +0.j        ])"
       ]
      }
     ],
     "prompt_number": 121
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}