Appendix B, Problem 16 - Calculate bubble point pressure of system with defined convergence pressure.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem 16\n",
    "Estimate bubble point pressure with modified Wilson correlation with 3 component mixture defined in Problem 15\n",
    "* Note the tabulated and reported results for this are confusing in the example as it appears that convergence pressure pk has been inadverdantly substituted for bubble point pressure. The actual problem ASSUMES pk = 2,000 psia (per Problem 15), and seeks to find the pressure at which Eq 3.165 is satisfied (being the bubble point pressure)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "from tabulate import tabulate\n",
    "R = 10.73146 # Universal gas constant, psia-ft3/R-lbm mol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Per Table B-18 and B-19\n",
    "comps = ['C1', 'C2', 'nC5']\n",
    "z = [0.20, 0.32, 0.48]\n",
    "pc = [667.8, 616.3, 488.6]\n",
    "tc = [343.0, 665.7, 845.4]\n",
    "w = [0.0115, 0.1454, 0.251]\n",
    "\n",
    "# As specified by problem scenario\n",
    "pk = 2000 # Convergence pressure psia \n",
    "p = 500 # psia\n",
    "t = 160 # deg F\n",
    "A2 = 0.7 # problem characterizes this as A0, but I assume this is a misprint as A0 doesn't factor into these equations\n",
    "pb = 2000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setup a function we can call for the modified Wilson correlation\n",
    "def k_mod_wilson(w, tc, pc, t, p, pk, A2):\n",
    "    A1 = 1-((p-14.7)/(pk-14.7))**A2 # Eq 3.160 - assuming pressures given in psia, Eq 3.160 is in terms of psig\n",
    "    tr = t/tc\n",
    "    pr = p/pc\n",
    "    return (pc/pk)**(A1-1)*math.exp(5.37*A1*(1+w)*(1-1/tr))/pr # Eq 3.159"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate all the Ki's\n",
    "ncomps = len(comps)\n",
    "K = [k_mod_wilson(w[i], tc[i], pc[i], t+460, p, pk, A2) for i in range(ncomps)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Fpbub(pbub, *args):\n",
    "    w, tc, pc, t, pk, A2 = args\n",
    "    ncomps = len(w)\n",
    "    K = [k_mod_wilson(w[i], tc[i], pc[i], t+460, pbub, pk, A2) for i in range(ncomps)]\n",
    "    return 1 - sum((z[i]*K[i]) for i in range(ncomps)) # Eq 3.165"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Table B-22 Pressure Guess Calculations (Problem 16)\n",
      "\n",
      "Bubble Point Pressure = 1375.1 psia\n",
      "\n",
      "Component      zi  Ki                   yi\n",
      "-----------  ----  ----------  -----------\n",
      "C1           0.2   1.9815      0.3963\n",
      "C2           0.32  0.9956      0.3186\n",
      "nC5          0.48  0.594       0.2851\n",
      "Total        1     1-Sum(yi):  2.22045e-16\n"
     ]
    }
   ],
   "source": [
    "from scipy.optimize import brentq\n",
    "# Calculate bubble point pressure consistent with specified convergence pressure pk\n",
    "data = (w, tc, pc, t, pk, A2)\n",
    "pbub = brentq(Fpbub,14.8,pk-1, args = data) # Use SciPy Python solver Brentq instead of Newton iterations, between standard conditions pressure and 1psi less than pk (Pbub = pk is trivial solution)\n",
    "\n",
    "K = [k_mod_wilson(w[i], tc[i], pc[i], t+460, pbub, pk, A2) for i in range(ncomps)]\n",
    "y = [K[i]*z[i] for i in range(ncomps)]\n",
    "\n",
    "# And print out summary\n",
    "print('Table B-22 Pressure Guess Calculations (Problem 16)\\n')\n",
    "print('Bubble Point Pressure =', round(pbub,1),'psia\\n')\n",
    "header=['Component', 'zi', 'Ki', 'yi',]\n",
    "table = []\n",
    "for i in range(ncomps):\n",
    "    table.append([comps[i], z[i], round(K[i],4), round(y[i],4)])\n",
    "table.append(['Total', sum(z), '1-Sum(yi):', 1-sum(y) ])\n",
    "print(tabulate(table,headers=header))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Note how Pbub has two solutions. The upper Pbub = pk is the trivial solution for when all Ki's = 1, hence why we set upper search range to 1psi less than pk*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEWCAYAAABBvWFzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd4VGXax/HvnR4IhBIJJfTewUSwsRIruCK6oqKuYl9dcYu6Lq51XV3rrq59sb1YEEGl2BZRiagovXdCTUJPIAQSQpL7/WMO65hNSJuZM5O5P9c1F6c8Z84vh8ncOe05oqoYY4wxdRXhdgBjjDH1gxUUY4wxPmEFxRhjjE9YQTHGGOMTVlCMMcb4hBUUY4wxPmEFxRhjjE9YQTGmCiKyRUQKRaTA6/VCLd/rHRHZISL5IrJeRG48Tts+IjJTRPaKSIU3jInIaBFZIyKHRCRTRIY40zuIiJbLfH9tMhtTXWI3NhpzfCKyBbhRVb/0wXv1Bjaq6hER6QFkAL9U1UUVtO0OnA7sBaapqpSbfw7wGnA5MB9oBaCq2SLSAdgMRKtqSV1zG1MdtodiTB2ISKSIPO3sRWwSkducPYOoitqr6ipVPXJs1Hl1rqTtOlV9HVhVyer/Cjysqj+qapmqZqtqdl1/JmNqywqKMXVzE3ABMBBIA0ZVtYCIvCQih4G1wA7gs5quVEQinfWdICIbRSRLRF4QkfhyTbc6894UkaSarseYmrCCYkz1TBOR/V6vm5zplwHPqup2Vc0FHqvqjVT1t0AjYAjwEXDk+EtUKBmIxlPAhgAD8BS1+5z5e4GTgPZAqrO+d2uxHmOqzQqKMdVzkao28Xq96kxvDWz3are1Om+mqqWq+h2QAtxaizyFzr/Pq+oOVd0L/BM433n/AlVdqKolqroLGAucKyKNa7EuY6rFCooxdbMDaOs13q6Gy0dRyTmU41HVPCALzzmYai3i/CvHbWVMHVhBMaZuJgO/E5EUEWkKjKusoYi0cC7zTXBO5p8HXAF8XUl7EZE4IMYZjxORWK8mbwK3O+/bFPgD8InTdrCIdBeRCBFpDjwHZKjqAR/8zMZUyAqKMdXzcbl7OqY6018FZgLLgMV4zolURvEc3soC8oCngT+o6nQAEWnnvPexvZz2eA5tHbvKqxBY5/V+fwMWAOuBNcAS4FFnXifgP8BBYCWe8zRX1OYHN6a67D4UY3zI7v8w4cz2UIwxxviEFRRjjDE+YYe8jDHG+ITtoRhjjPGJCvsbqq+SkpK0Q4cOtVr20KFDNGzY0LeB/CSUskJo5Q2lrBBaeUMpK4RW3rpkTUpKYubMmTNVdViVjVU1bF6pqalaW7Nnz671soEWSllVQytvKGVVDa28oZRVNbTy1jUrsFCr8R1rh7yMMcb4hBUUY4wxPmEFxRhjjE9YQTHGGOMTVlCMMcb4hBUUY4wxPmEFxRhjjE+E1Y2NxgSzw8UlZOcVsufgEQ4Vl1J4tJTC4hIOF5dyuLiU4pIyYqMjiI+O9LxiIomLjqRBTCQtGsXRukkcjeKi3f4xTBhztaCIyDDgX0Ak8JqqPl5ufizwFp5nYu8DLlfVLV7z2wGrgYdU9elA5Tamto6WlrF+10GWZx0gc3cB2fsLycorJHt/IbmHiuv8/o3iomjTJJ7WTeJp3SQOPXCUhC259GzVmIax9vej8S/XPmEiEgm8CJyD54FDC0Rkhqqu9mp2A5Cnql1EZDTwBHC51/xngM8DldmYmtq67xCLtuaxPOsAy7L2szonnyMlZQDERkXQpmk8KU0b0DclkTZN4klpGk+LRnE0jPXsecTHRNHA2RuJiYyguLSMQmfvpeio59/DxaXsyi8iO6+QnP2FZO8vImd/IYu25nGg8CjvrvkBEeiU1JA+bRLp0zqRge2a0L9tE6Ij7ai38R03/2QZBGxU1U0AIjIJGIlnj+OYkcBDzvAHwAsiIqqqInIRsAk4FLjIxhxf0dFS5m/O5eu1u8lYt5st+w4DEB8dSZ82jfn1ye3pl5JI/5QmtGvWgIiImj3iPS7Cc5iraTXaqirTZs6mUbverMw5wMrsfOZvzmX60hwAGsVGcUrn5gzpdgK/6JpE++ah0S+VCV6udV8vIqOAYap6ozN+NTBYVcd6tVnptMlyxjOBwXgehfolnr2bu4CCyg55icjNwM0AycnJqZMmTapV3oKCAhISEmq1bKCFUlYIrbwVZS0oVhbsLGHZnlJW55ZSXArREdCzeST9kiLp3iyS1g2FyBoWD3/lzT+irMsrZdXeUlbsLWVfkec74IR4oe8JkQxuGUXXphFESGDzhtLnAEIrb12zpqenL1LVtKraubmHUtGntXx1q6zNX4FnVLVAqvjQq+p4YDxAWlqaDh06tOZJgYyMDGq7bKCFUlYIrbzHspaUljFnwx4+WJTFl6t3U1xaRkrTeEYPakF69xac0rk5cdGRbsetctuqKlv2HebbDXuYs34v32/cy9fbimiVGMeI/q25sH9rerduTFW/Z4HIGmxCKW+gsrpZULKAtl7jKUBOJW2yRCQKSARy8eyljBKRJ4EmQJmIFKnqC/6PbcJZdkEZj322ho+WZLPn4BGaNYzh1ye355LUNvRqFZgvXl8SETomNaRjUkOuOaUDh46U8OWaXcxYmsMb321m/JxNdDqhISP7t+Hyk9rSMjHO7cgmiLlZUBYAXUWkI5ANjAauLNdmBjAG+AEYBXztdKU85FgDEXkIzyEvKybGL1SVORv28u9vMpmbWUhUxGbSe7RgVGoK6d1bEBNVf05sN4yNYuSANowc0Ia8Q8V8vnInM5Zl8+xX63nu6w0M692SMad24KQOTUOueBr/c62gqGqJiIwFZuK5bPgNVV0lIg/j6Xt/BvA68LaIbMSzZzLarbwm/JSUlvHpih288s0m1uzIp2XjOC7rHs3dl55BUkKs2/H8rmnDGK4c3I4rB7dje+5h3v5xK+8v2M6nK3bQs1VjxpzSnpED2hAf4/6hPRMcXL0wXVU/Az4rN+0Br+Ei4NIq3uMhv4QzYetwcQnvL9jOa99uJnt/IV1aJPDUqH6MHNCGud/NCYtiUl7bZg34y/k9+ePZ3Zi2NJsJc7cw7qMVPPb5Wq49tQPXn96RxHi7qTLc2Z1OxjhKy5QpC7fzz1nr2X3wCGntm/LXC3tzZo8WNb68t76Kj4nkikHtGH1SW+ZvzuX17zbzr6828MZ3m7nu9I7ccFpHEhtYYQlXVlBM2FNVMtbt4bHP17B+VwEntmvC81cMZHCn5m5HC1oiwuBOzRncqTmrc/J57qsNPPfVBt60whLWrKCYsLYy+wB//2wNczP30b55A1666kSG92lpJ5xroFfrxrxyder/FJabftGJm4Z0snMsYcQKiglLuYeKefTTNXy4OIumDaJ5cEQvrhrcvl5dsRVoxwrLmh35PDNrPf+ctZ5J87cx7vyejOjXyop0GLCCYsKKqjJjWQ5//Xg1+YVH+c0ZnbgtvQuNrZden+nZqjHjr0njx037+OvHq/nde0t4a+4WHhjRi34pTdyOZ/zICooJG1l5h7lv2koy1u1hQNsmPH5JX3q0bOx2rHrr5E7N+eT205mycDtPf7GOC1/4nlGpKdx9XndaNLYbJOsjKyim3istUybM3cLTX6wD4MERvbjmlA6u9K0VbiIjhNGD2vHLfq14YfZG3vxuCzNX7mTc+T244qR2dvVcPWMFxdRr2/Yd5vfvL2HJtv0M7X4Cj1zUh5SmDdyOFXYaxUVzz/CeXHFSO/4ydQX3Tl3JjKU5PH5JPzomWS/H9YWdgTT11vSl2Zz/3Ldk7i7g2csH8Oa1J1kxcVmHpIa8e+NgnrikL6t35DPs2Tm88k0mJaVlbkczPmB7KKbeOVxcwkMzVjF5YRap7Zvyr9EDrJAEERHh8pPaMbR7Cx6YvpLHP1/LJ8tzuLR9qdvRTB1ZQTH1yuqcfG5/bzGb9h5ibHoX/nB2V6LsqYRBKblxHP++Oo3PV+zg/umrePiHIxxqvJHf/KKznd8KUfabZuoFVeWtH7Zw0Uvfc7CohHdvGMxd53W3YhIChvdtxVd3nEFqciRP/mcdv35tHjsPFLkdy9SC/baZkFd0tJQ7Jy/jgemrOL1LEp//fgindklyO5apgcQG0dzaP5YnR/VjWdZ+hv1rDv9ZudPtWKaGrKCYkLY7v4jLx//IR0uyufOcbrw+Jo3mYdgbcH0gIlyW1pZPfzeEds0acMs7i/jL1BUUFtu5lVBhBcWErOVZ+xnxwnds2HWQf1+dyu1ndbXuPeqBjkkN+eCWU7nljM68N38bFzz/Let3HXQ7lqkGKygmJE1fms2lr/xAVEQEH956Kuf1bul2JONDMVERjBveg3dvGEx+UQkXvfg9nywv/4RwE2ysoJiQUlamPDVzLb+ftJT+KU2YMfY0eray7lPqq1O7JPHJ7afTs1Vjxk5cwt8/W2P3rAQxKygmZBSXlPG7SUt4cXYmo09qyzs3DrbzJWEguXEc7910Mlef3J7xczZxzRvz2VdwxO1YpgJWUExIOHSkhBsmLOCT5TsYN7wHj/2qr3U1H0ZioiL420V9eGpUPxZuzWPE89+xPGu/27FMOfYbaYJe3qFirnptHt9v3MuTo/pxyxmd7eR7mLo0rS0f3nIqIsKoV37go8VZbkcyXqygmKC240Ahl/77B1bvyOflX6dyWVpbtyMZl/VNSeTj208ntV1T7pi8jGdmrUdV3Y5lcLmgiMgwEVknIhtFZFwF82NF5H1n/jwR6eBMHyQiS53XMhG5ONDZjf9l7ilg1Ms/sPNAEROuG2RXcpn/atYwhgnXD2JUagr/+moDd0xexpESu1/Fba715SUikcCLwDlAFrBARGao6mqvZjcAearaRURGA08AlwMrgTRVLRGRVsAyEflYVUsC/GMYP1mRdYAxb84nQmDSzSfTp02i25FMkImJiuCpUf1o36wB/5i1nuz9hYy/OpUmDWLcjha23NxDGQRsVNVNqloMTAJGlmszEpjgDH8AnCUioqqHvYpHHGD7u/XIiqwDXPXaj8RHRzLlllOtmJhKiQi3n9WVf40ewNJt+/nVS3PZsveQ27HClrh17FFERgHDVPVGZ/xqYLCqjvVqs9Jpk+WMZzpt9orIYOANoD1wtapOrWQ9NwM3AyQnJ6dOmjSpVnkLCgpISEio1bKBFkpZ4ed5t+aX8uSCIuIihXsGx5EUH1yn+UJ52wa7umZdn1fKvxYXEQH87sQ4ujaN9F24CoTTtk1PT1+kqmlVNlRVV17ApcBrXuNXA8+Xa7MKSPEazwSal2vTE5gPxFW1ztTUVK2t2bNn13rZQAulrKo/5V2dc0AH/HWmnvL3L3XbvkPuhqpEqG7bUOCLrJv2FOgZT36t3e/7TGev3VX3UMcRTtsWWKjV+F5388+/LMD7kp0UoHzfCv9tIyJRQCKQ691AVdcAh4A+fktq/G7dzoNc9do84qIjee/mk2nbzB6IZWquY1JDPrj1VDolJXDTWwv5eJl11xJIbhaUBUBXEekoIjHAaGBGuTYzgDHO8Cjga1VVZ5koABFpD3QHtgQmtvG17IIyrnz1R6IjhYk3nUz75vaMcVN7SQmxTPrNyQxo24TfTVrCxHnb3I4UNlwrKOo5qT4WmAmsASar6ioReVhELnSavQ40F5GNwB3AsUuLT8dzZddSYCrwW1XdG9ifwPjCxt0FPDG/iIgITzHpmGTFxNRd47ho3rp+MGd0O4G/TF3ByxmZbkcKC64+AlhVPwM+KzftAa/hIjznWsov9zbwtt8DGr/K3l/Ir1+bByjv3TSYzieExglOExriYyIZf3Uad05ZxhP/WcuBwqP8eVh362XBj+yZ8sYV+w8XM+aN+RwqLuHuk+Lp0qKR25FMPRQTFcGzlw+gcVwUr3yTSX7RUR4Z2YcIe2a9X1hBMQFXdLSUGyYsZNu+w7x1wyCKtq1wO5KpxyIjhEcu6kPj+GhezshEVXn0or5WVPzACooJqJLSMsZOXMLibXm8eOWJnNypORl2ztT4mYhw93ndiRThhdkbAayo+IEVFBMwqsr901fx5ZpdPDSiF+f3beV2JBNGRIQ7z+0GYEXFT6ygmID511cbeG/+Nn47tDPXntbR7TgmDB0rKory4uxMQHj0Ijun4itWUExATJy3jWe/3MAlJ6bwp/O6ux3HhDER4a5zPZ9BT1HBioqPWEExfjdn/R7um7aCM7qdwOOX9LXLNo3rrKj4hxUU41eb9hQwduJiuiU34sWrTiQ6Mrg6ezTh61hRUYWXMjJpEBPJfb/saX/w1IEVFOM3BwqPcuOEhURFRvDqNWkkxNrHzQQXEeFP53XncHEpr3+3mcT4aH53Vle3Y4Us+w03flFaptz+3hK25R7m3RsHW2ePJmiJCA9c0IuDRSX8c9Z6GsdF2UUjtWQFxfjFY5+tYc76PTz2q74M7tTc7TjGHFdEhPDEJX05WHSUhz5eTeP4aH51YorbsUKOHdA2Pjdl4XZe+24zY05pzxWD2rkdx5hqiYqM4LkrBnJq5+b86YPlfLFqp9uRQo4VFONTi7bmcu/UlZzWpTn3X9DL7TjG1EhcdCTjr0mjT5tExr63hLmZ1ol5TVhBMT6z40Ahv3l7Ea2bxPHilScSZVd0mRCUEBvFhOtOokPzBtw0YSErsg64HSlk2G+88YmjTh9dhcWlvDYmjSYNYtyOZEytNWkQw9s3DKZpwxiu+78FbM897HakkGAFxfjEk/9Zy6KteTx+ST/rit7UC8mN4/i/607iaGkZY96cz/7DxW5HCnpWUEydfbFqJ69+u5lrTmnPiP6t3Y5jjM90adGIV69JIyu3kJveWkjR0VK3IwU1KyimTrbtO8ydU5bRLyWRe3/Z0+04xvjcoI7N+Mdl/VmwJY87pyyjrEzdjhS07D4UU2tHSkq5beJiBHjxyhOJjYp0O5IxfjGif2ty9hfy2OdrSWkSzz3n2x9PFbGCYmrtkU/WsCL7AK9ek2Z3wpt67+ZfdCJ7fyH/nrOJNk3jsTus/perh7xEZJiIrBORjSIyroL5sSLyvjN/noh0cKafIyKLRGSF8++Zgc4e7mYsy+HtH7dy8y86cU6vZLfjGON3IsKDI3pzds9kHpqxiiW7S9yOFHRcKygiEgm8CAwHegFXiEj5O+FuAPJUtQvwDPCEM30vMEJV+wJjgLcDk9oAZO4p4J4Pl5PWvqk928SElcgI4fkrBtK3TSKvLDvCmh35bkcKKm7uoQwCNqrqJlUtBiYBI8u1GQlMcIY/AM4SEVHVJaqa40xfBcSJSGxAUoe54pIyfj9pCTFRETx/5UDrjt6EnfgYz9308VHCjRMWsq/giNuRgoaounPFgoiMAoap6o3O+NXAYFUd69VmpdMmyxnPdNrsLfc+t6jq2ZWs52bgZoDk5OTUSZMm1SpvQUEBCQkJtVo20PyZ9YP1xXyy6Si3D4wlNdk3p+Bs2/pPKOUNpawAq3YU8OwKoWNiBHefFEdUED+cq67bNj09fZGqplXZUFVdeQGXAq95jV8NPF+uzSogxWs8E2juNd7bmda5OutMTU3V2po9e3atlw00f2Wdt2mfdhj3if75g2U+fV/btv4TSnlDKauqJ++Mpdna/s+f6J+mLNWysjK3I1WqrtsWWKjV+I5183hFFtDWazwFyKmsjYhEAYlArjOeAkwFrlHVTL+nDXP5RUf54/tLadesgXX6aIxjRP/W/O7MLkxemMUb329xO47r3CwoC4CuItJRRGKA0cCMcm1m4DnpDjAK+FpVVUSaAJ8C96jq9wFLHMYemr6KnflFPHP5ABrakxeN+a8/nN2NYb1b8uinq/lm/R6347jKtYKiqiXAWGAmsAaYrKqrRORhEbnQafY60FxENgJ3AMcuLR4LdAHuF5GlzqtFgH+EsPHxshw+WpLN7Wd24cR2Td2OY0xQiYgQ/nl5f7q3bMzYiYvZuLvA7UiucfUSHVX9TFW7qWpnVX3UmfaAqs5whotU9VJV7aKqg1R1kzP9EVVtqKoDvF673fxZ6qsdBwq5d+oKBrRtwtj0Lm7HMSYoNYiJ4rUxacRGRXDTWwvJLzrqdiRX2DWfplJlZcqdk5dRUqY8e/kAe76JMcfRpkk8L12Vyvbcw9w1OTz7/LJvCFOpN77fzNzMfTw4ohcdkhq6HceYoDeoYzP+cn5Pvli9i1fmhN+1QlZQTIU27SngqZnrOLtnMpelta16AWMMANed1oER/Vvz9Mx1fLchvB4hbAXF/I+yMuXPHy4nNiqCv1/cB5HgvWHLmGAjIjxxSV+6tmjE7e8tJisvfJ72aAXF/I935m1lwZY87r+gFy0ax7kdx5iQ0yAmileuTqWkVPntu4vD5sFcVlDMz2zPPczjn6/lF91OYFRqittxjAlZHZMa8o/L+rM86wAPzVjldpyAsIJi/ktV+cvUFQjYoS5jfODc3i25Lb0zkxZsZ9L8bW7H8TsrKOa/pizK4tsNexk3vAcpTe2BWcb4wh3ndGdI1yQemLGKVTkH3I7jV1ZQDAC78ot45JPVDOrYjKsGt3c7jjH1RmSE8OzlA2jaIJqxE5dQcKT+PpjLCopBVblv2kqOlJTxxCX9iAjibriNCUXNE2J5bvRAtu47xL1TVxzrLb3esYJi+GT5Dmat3sWd53ajo93AaIxfDO7UnD+e3Y3pS3OYvHC723H8wgpKmNt/uJiHZqyif0oi15/W0e04xtRrv03vwmldmvPgjFWs23nQ7Tg+ZwUlzD01cx15h4t57Ff9rK8uY/wsMkJ45vIBJMRGc9vExRwurl/nU+wbJIwt3b6fifO3ce2pHenVurHbcYwJCy0axfHs5QPI3FPAg9Pr1/0pVlDCVGmZcu/UFbRoFMsfz+nqdhxjwsrpXZMYm96FKYuy+GhxlttxfKbaj94TkTRgCNAaKARWAl+qaq6fshk/eufHrazKyeeFKwfSKC7a7TjGhJ3fn9WVeZtzuW/aSga0bUKnExLcjlRnVe6hiMi1IrIYuAeIB9YBu4HTgVkiMkFE2vk3pvGl3QeLeHrmOoZ0TeKXfVu5HceYsBQVGcFzowcSExXBH95fytHSMrcj1Vl19lAaAqepamFFM0VkANAVqP/9CtQTf/90DUdKyvjrhb2texVjXNQyMY7HLu7Lre8u5tkv1/On83q4HalOqtxDUdUXKysmzvylqvqVb2MZf5mbuZdpS3O45YxO9WIX25hQN7xvKy5LS+GljEzmbdrndpw6qfZJeRHpJCIfi8heEdktItNFpJM/wxnfKi4p4/5pK2nbLJ7f2vPhjQkaD47oTftmDbhj8jIOFIbu8+hrcpXXRGAy0BLPifkpwHv+CGX847XvNpG55xAPX9iHuOhIt+MYYxwNY6N45vIB7Mwv4oHpK92OU2s1KSiiqm+raonzegeoU4c0IjJMRNaJyEYRGVfB/FgRed+ZP09EOjjTm4vIbBEpEJEX6pIhXGTvL+S5rzZwXu9k0nu0cDuOMaacge2a8oezujJ9aQ7TlmS7HadWqnOVVzMRaQbMFpFxItJBRNqLyN3Ap7VdsYhEAi8Cw4FewBUi0qtcsxuAPFXtAjwDPOFMLwLuB+6q7frDzROfr0UV7r+g/CY2xgSL36Z3Ia19U+6ftpLtuaH36ODq7KEsAhYClwO/AWYDGcCtwHV1WPcgYKOqblLVYmASMLJcm5HABGf4A+AsERFVPaSq3+EpLKYKi7bmMmNZDr/5RSd7zokxQexY1ywAd0xeSmlZaPVKLG51oywio4BhqnqjM341MFhVx3q1Wem0yXLGM502e53xa4E072UqWM/NwM0AycnJqZMmTapV3oKCAhISQuOqKO+sZao88mMRuUXK40PiiYsKvsuEQ3XbhoJQyhtKWcG/eefmlDB++RFGdY3mgs4xdX6/umZNT09fpKppVbWryZ3y11Q0XVXfqkkw77es6O1q0ea4VHU8MB4gLS1Nhw4dWpPF/ysjI4PaLhto3lk/WpzFpgPL+Mel/RkWpM+ID9VtGwpCKW8oZQX/5j1DlWxdwvTVO7nxlyfTo2Xd+toL1LatyUn5k7xeQ4CHgAvrsO4soK3XeAqQU1kbEYkCEgHr6qWaDheX8MR/1tIvJZGLB7ZxO44xpppEhL9d1IfE+GjunLwsZO6ir3ZBUdXbvV43AQOBuuyLLQC6ikhHEYkBRgMzyrWZAYxxhkcBX2t9fdSZH7zyzSZ25R/hgQt62VMYjQkxzRrG8OjFfT197n290e041VKX3oYP4+lypVZUtQQYC8wE1gCTVXWViDwsIsf2fF4HmovIRuAO4L+XFovIFuCfwLUiklXBFWJhLXt/If/+JpMR/VuT1qGZ23GMMbVwXu+WXDywDS/O3siKrANux6lSTc6hfMxP5y8igZ54bnSsNVX9DPis3LQHvIaLgEsrWbZDXdZd3z35n7UA/HlYd5eTGGPq4qERvZmbuZc7pyzl49tPJzYqeG9KrnZBAZ72Gi4Bth67+soEl415pUxfmsPtZ3axy4SNCXGJDaJ5/JJ+XPfmAp6ZtYFxw4O3A8manEP5Bk/X9YlAMzxFxQSZsjJl4tpiWjSK5ZYzOrsdxxjjA+ndWzD6pLaMn5PJ4m15bsepVE06h7wRmA/8Cs8J8h9F5Hp/BTO18/HyHDYdKOPuYT1oGFuTHVBjTDC795c9aZUYz12Tl1FYXOp2nArV5KT8n4CBqnqtqo4BUoE/+yeWqY3ikjKe/mIdbRtF8Cu7TNiYeqVRXDRPjerHpr2HeGrmOrfjVKgmBSULOOg1fhDY7ts4pi4mztvK9txCLu0WbZcJG1MPndoliV+f3I43524OykNf1ekc8g4RuQPIBuaJyEMi8iDwIxAaF0eHgYIjJTz/9UZO7tSMvknBexWIMaZu/jysBy0bxzHuw+UUlwTXDY/V2UNp5LwygWn8dOnwdGCHn3KZGnrt203sO1TMuOE97bG+xtRjjeKiefTiPqzfVcBLGcH1N32VZ21V9a8icgLQHk/vwPv9H8vUxN6CI7w6ZxPD+7RkQNsmZGS6ncgY409n9kjmogGteXH2Rob3aUX3lo3cjgRU75DXDcAq4Hlgrddd7CZIvPD1RopKyrjrPLttZt4oAAAU70lEQVSJ0Zhw8cCI3jSKi+bPHy4Pmm7uq3PI649Ab1U9BTgVuMe/kUxNbM89zLvztnJZWgqdTwidrr+NMXXTrGEMD47oxdLt+/m/uVvcjgNUr6AUq+oeAFXdBMT6N5KpiX/OWk+ECL8/q5vbUYwxAXZh/9ac1aMFT89cx7Z97j/hsToFJUVEnjv2qmDcuGTNjnymLc3mutM60jIxzu04xpgAExEeubgPkRHCX6auwO3O2KtzK/Wfyo0v8kcQU3NP/mctjeOiudW6WDEmbLVKjGfc8B7cN20lUxZmcdlJbateyE+qc5XXhKramMCbt2kfs9ft4Z7hPUhsEO12HGOMi64c1I4Zy3J49LM1nNmzBUkJ7pyZqM5VXuNFpE8l8xqKyPUicpXvo5nKqCpPf7GOlo3jGHNqB7fjGGNcFhEh/P3ivhwuLuHRT9e4l6MabV4CHhCRNSIyRUReEpE3RORbYC6emx4/8GtK8zPfb9zHgi153JbembhouyveGANdWiRwyxmdmbokm+837nUlQ3UOeS0FLhORBCANaAUUAmtUNTh7KKvHVJVnvlxP68Q4V4+VGmOCz23pXZixLIf7pq3k898PCfgfnNU55NUOQFULVDVDVd9T1WlWTNzx7Ya9LNqax2/TuwT1k9uMMYEXFx3JIxf1YfPeQ7zsQpcZ1TnkNe3YgIh86McspgrH9k7aNInnsjTbOzHG/K8hXU/gwv6teTkjk8w9BQFdd3UKindPg538FcRU7Zv1e1iybT+3pXchJqomTx4wxoST+y7oSWx0BPdPWxnQe1Oq862klQzXmYgME5F1IrJRRMZVMD9WRN535s8TkQ5e8+5xpq8TkfN8mSsYefZONtCmSTyjUlPcjmOMCWItGsVx97AezM3cx9Ql2QFbb3UKSn8RyReRg0A/ZzhfRA6KSH5tVywikcCLwHCgF3CFiPQq1+wGIE9VuwDPAE84y/YCRgO9gWHAS8771VsZ6/awbPt+bj/T9k6MMVW7alA7BrRtwqOfrqGgODB7KVV+M6lqpKo2VtVGqhrlDB8bb1yHdQ/C0x3+JlUtBiYBI8u1GQkcu7HyA+As8TzsYyQwSVWPqOpmPA/6GlSHLEHt2LmTts3iucT2Towx1XDs3pT9hUeZvL44IOusTtcr/tKGnz9COAsYXFkbVS0RkQNAc2f6j+WWrfAh6iJyM3AzQHJyMhkZGbUKW1BQUOtl62rJ7hKWZx3h+j4xfP/tnCrbu5m1NkIpbyhlhdDKG0pZIXTyntMukllbj/Lh51/TPN6/RzfcLCgVPVaw/H5ZZW2qs6xnoup4YDxAWlqaDh06tAYRf5KRkUFtl60LVeXpF76jXbNI7rniDKIjq/5AuJW1tkIpbyhlhdDKG0pZIXTynnRKCVO/mMMlw8/0+7rcPBifBXhf+5oC5FTWRkSigEQgt5rL1guzVu9iZXY+t5/ZpVrFxBhjvDWMjSKlUWC+O9z8hloAdBWRjiISg+ck+4xybWYAY5zhUcDX6rkGbgYw2rkKrCPQFZgfoNwBo6o89/UG2jdvwMUDKzyiZ4wxQcO1Q17OOZGxwEwgEnhDVVeJyMPAQlWdAbwOvC0iG/HsmYx2ll0lIpOB1UAJcJuqlrryg/jRnA17WZmdz5OX9CPK9k6MMUHOzXMoqOpnwGflpj3gNVwEXFrJso8Cj/o1oMtenL2RVolxXGR7J8aYEGB/9gaphVtymb85l5uGdLL7TowxIcG+qYLUSxmZNG0QzehB1meXMSY0WEEJQqtz8vl67W6uP60jDWJcPSppjDHVZgUlCL38TSYJsVFcc0oHt6MYY0y1WUEJMpv3HuLT5TlcdXI7e1a8MSakWEEJMv/+JpOoyAhuOL2j21GMMaZGrKAEkZ0HivhwcRaXpaXQolGc23GMMaZGrKAEkVe/3USZwm9+0dntKMYYU2NWUIJE3qFiJs7bxoX9W9O2WQO34xhjTI1ZQQkSb87dQuHRUm4dansnxpjQZAUlCBw6UsKEuVs4t1cy3ZIbuR3HGGNqxQpKEJiycDsHCo/ymzM6uR3FGGNqzQqKy0rLlDe+38LAdk1Ibd/M7TjGGFNrVlBcNmv1TrblHuamIbZ3YowJbVZQXPbqt5tp2yye83q3dDuKMcbUiRUUFy3elseirXlcf1pHIiPE7TjGGFMnVlBc9Pq3m2kcF8VladZFvTEm9FlBccn23MN8vnIHVw5uT8NY66LeGBP6rKC45I3vNxMhwrWndnA7ijHG+IQVFBccKDzK5AXbubB/a1omWieQxpj6wQqKCybN38ah4lJuGGJd1Btj6g9XCoqINBORWSKywfm3aSXtxjhtNojIGK/pj4rIdhEpCFxq3zhaWsb/zd3CqZ2b07t1ottxjDHGZ9zaQxkHfKWqXYGvnPGfEZFmwIPAYGAQ8KBX4fnYmRZyPl2+gx0HiuxGRmNMveNWQRkJTHCGJwAXVdDmPGCWquaqah4wCxgGoKo/quqOgCT1IVXl1W830aVFAmd0O8HtOMYY41OiqoFfqch+VW3iNZ6nqk3LtbkLiFPVR5zx+4FCVX3aq02BqiZUsa6bgZsBkpOTUydNmlSrzAUFBSQkHHdVVVqbW8rj84u4rncMZ7T13/PifZE1kEIpbyhlhdDKG0pZIbTy1jVrenr6IlVNq6qd326AEJEvgYr6E7m3um9RwbQaVz9VHQ+MB0hLS9OhQ4fW9C0AyMjIoLbLHjP53UU0abCPP48+k7joyDq91/H4ImsghVLeUMoKoZU3lLJCaOUNVFa/FRRVPbuyeSKyS0RaqeoOEWkF7K6gWRYw1Gs8BcjwacgAytlfyMxVu7hxSEe/FhNjjHGLW+dQZgDHrtoaA0yvoM1M4FwRaeqcjD/XmRaS3p23FVXl6pPbux3FGGP8wq2C8jhwjohsAM5xxhGRNBF5DUBVc4G/AQuc18PONETkSRHJAhqISJaIPOTCz1BtRUdLeW/+ds7umUxKU3tevDGmfnKlEylV3QecVcH0hcCNXuNvAG9U0O5u4G5/ZvSlT5fvIPdQMWOsmxVjTD1md8r7maoy4YctdGmRwKmdm7sdxxhj/MYKip8t3b6f5VkHGHNKe0TsmSfGmPrLCoqfTZi7hUaxUfzqxBS3oxhjjF9ZQfGj3QeL+HTFDkalpdgzT4wx9Z4VFD+aNH87R0vtUmFjTHiwguInR0vLeHfeVs7odgKdTgiN7hmMMaYurKD4yX9W7mRX/hHGnGp7J8aY8GAFxU/e+mEL7Zo1YGi3Fm5HMcaYgLCC4gercg6wYEse15zSnogIu1TYGBMerKD4wTs/biMuOoJLU9u6HcUYYwLGCoqPFRwpYcbSbEb0a01iA/8988QYY4KNFRQfm740m0PFpVwxuJ3bUYwxJqCsoPiQqjJx3jZ6tGzEwLZNql7AGGPqESsoPrQ86wCrcvK5anA767fLGBN2rKD40MR524iPjmTkwDZuRzHGmICzguIj+UVHmbEshwv7t6ZxnJ2MN8aEHysoPjJ9STaFR0u50k7GG2PClBUUH1BV3p23jd6tG9MvJdHtOMYY4worKD6wZPt+1u48yJV2Mt4YE8asoPjAxHnbaBgTycgBdjLeGBO+XCkoItJMRGaJyAbn36aVtBvjtNkgImOcaQ1E5FMRWSsiq0Tk8cCm/7kDhUf5ZHkOFw5oQ4I9RMsYE8bc2kMZB3ylql2Br5zxnxGRZsCDwGBgEPCgV+F5WlV7AAOB00RkeGBi/6+pi7MoOlrGVXYy3hgT5twqKCOBCc7wBOCiCtqcB8xS1VxVzQNmAcNU9bCqzgZQ1WJgMeDKA9tVlYnzt9G3TSJ92tjJeGNMeBNVDfxKRfarahOv8TxVbVquzV1AnKo+4ozfDxSq6tNebZrgKShnq+qmStZ1M3AzQHJycuqkSZNqlbmgoICEhJ8/eXFDXimPzivi2t4xDG0bPPeeVJQ1mIVS3lDKCqGVN5SyQmjlrWvW9PT0RaqaVlU7vx30F5EvgZYVzLq3um9RwbT/Vj8RiQLeA56rrJgAqOp4YDxAWlqaDh06tJqr/7mMjAzKL/vx5GUkxO7kT5el0zCIzp9UlDWYhVLeUMoKoZU3lLJCaOUNVFa/fQuq6tmVzRORXSLSSlV3iEgrYHcFzbKAoV7jKUCG1/h4YIOqPuuDuDV2sOgon63YwUUDWwdVMTHGGLe4dQ5lBjDGGR4DTK+gzUzgXBFp6pyMP9eZhog8AiQCfwhA1gp9unwHhUdLuTTNHqJljDHgXkF5HDhHRDYA5zjjiEiaiLwGoKq5wN+ABc7rYVXNFZEUPIfNegGLRWSpiNwY6B9gyqIsurRIsG7qjTHG4cqxGlXdB5xVwfSFwI1e428Ab5Rrk0XF51cCZuPuAhZtzeOe4T3sznhjjHHYnfK1MGXRdiIjhItPtDvjjTHmGCsoNVRSWsZHi7NJ796CFo3i3I5jjDFBwwpKDX2zfg97Dh7h0jRX7qU0xpigZQWlhqYszCIpIYYze7RwO4oxxgQVKyg1sK/gCF+u2cXFA9sQHWmbzhhjvNm3Yg1MXZJNSZnavSfGGFMBKyjVpKp8sCiL/m2b0C25kdtxjDEm6FhBqaYt+WWs3XmQy+xkvDHGVMgKSjV9m1VCbFQEI/q3djuKMcYEJSso1VB0tJQfd5QwvE9LGscFTzf1xhgTTKygVMPMVTs5XAKX2cl4Y4yplBWUapiyMIukeOHkTs3djmKMMUHLHuRRBVWlR8tGtI3KJyLCOoI0xpjK2B5KFUSE+y7oxXkd7NyJMcYcjxUUY4wxPmEFxRhjjE9YQTHGGOMTVlCMMcb4hBUUY4wxPmEFxRhjjE9YQTHGGOMTVlCMMcb4hKiq2xkCRkT2AFtruXgSsNeHcfwplLJCaOUNpawQWnlDKSuEVt66ZN0LoKrDqmoYVgWlLkRkoaqmuZ2jOkIpK4RW3lDKCqGVN5SyQmjlDVRWO+RljDHGJ6ygGGOM8QkrKNU33u0ANRBKWSG08oZSVgitvKGUFUIrb0Cy2jkUY4wxPmF7KMYYY3zCCooxxhifCNuCIiJviMhuEVnpNe1vIrJcRJaKyBci0tqZLiLynIhsdOaf6LXMGBHZ4LzGBDKv17y7RERFJCkY8laybR8SkWxn2y4VkfO95t3jZF0nIud5TR/mTNsoIuP8kbWyvM702531rxKRJ4MhbyXb9n2v7bpFRJYGQ9bj5B0gIj86eReKyCBnejB+bvuLyA8iskJEPhaRxl7z3PwctBWR2SKyxvl8/t6Z3kxEZjnbaZaINHWmB2bbqmpYvoBfACcCK72mNfYa/h3wijN8PvA5IMDJwDxnejNgk/NvU2e4aaDyOtPbAjPx3LCZFAx5K9m2DwF3VdC2F7AMiAU6AplApPPKBDoBMU6bXgH8LKQDXwKxzniLYMhb2efAa/4/gAeCIetxtu0XwHCvz2pGEH9uFwBnOMPXA38Lhm0LtAJOdIYbAeudTE8C45zp44AnArltw3YPRVXnALnlpuV7jTYEjl2xMBJ4Sz1+BJqISCvgPGCWquaqah4wC6jyblJf5XU8A9ztldX1vMfJWpGRwCRVPaKqm4GNwCDntVFVN6lqMTDJaetzleS9FXhcVY84bXYHQ97jbVsREeAy4L1gyHqcvAoc+0s/Ecjxyhtsn9vuwBxneBZwiVdWNz8HO1R1sTN8EFgDtHHWNcFpNgG4yCuv37dt2BaUyojIoyKyHbgKeMCZ3AbY7tUsy5lW2fSAEJELgWxVXVZuVlDmBcY6u9tvHNsVP04mt7N2A4aIyDwR+UZETnKmB2tegCHALlXd4IwHa9Y/AE85v2dPA/c404Mx70rgQmf4UjxHBDhOpoBnFZEOwEBgHpCsqjvAU3SAFoHMawWlHFW9V1XbAu8CY53JUlHT40z3OxFpANzLT0XvZ7MrmOZqXuBloDMwANiB59AMx8nkZlaAKDyHAE4G/gRMdvYAgjUvwBX8tHcCwZv1VuCPzu/ZH4HXnenBmPd64DYRWYTn0FKxMz0osopIAvAh8IdyR1j+p2kF03ye1wpK5Sby0+5tFj/9ZQKQgmc3vbLpgdAZz7HbZSKyxVn3YhFpeZxcruVV1V2qWqqqZcCreA4NEIxZvXJ95BwimA+U4elgLyjzikgU8Cvgfa/JQZkVGAN85AxPIYg/C6q6VlXPVdVUPMU6M1iyikg0nmLyrqoe2567nENZOP8eO1QbmLy+PlkUSi+gAz8/AdfVa/h24ANn+Jf8/ITWfP3phNZmPH/JNnWGmwUqb7l5W/jppLzreSvYtq28hv+I5/gzQG9+fnJzE54Tm1HOcEd+OrnZO4CfhVuAh53hbngOC0gw5K3oc4DnuPc35aa5nrWSbbsGGOoMnwUsCuLP7bGLMSKAt4Drg2HbOtvoLeDZctOf4ucn5Z8M5Lb1ywcoFF54/trYARzFU6VvwFPtVwLLgY+BNl7/eS/i+etkBZDm9T7X4zkhtxG4LpB5y83fwk8FxdW8lWzbt50sy4EZ/LzA3OtkXYdz9Y8z/Xw8V69kAvcG+LMQA7zjfB4WA2cGQ97KPgfA/wG3VNA+GLft6cAiPF+284DUIP7c/t7ZTuuBx3F6F3F72zrbUJ3fp6XO63ygOfAVsMH5t1kgt611vWKMMcYn7ByKMcYYn7CCYowxxiesoBhjjPEJKyjGGGN8wgqKMcYYn7CCYoyfiEip06PuShGZIiINRKRD+V6Nq/E+D4nIXf7KaYyvWEExxn8KVXWAqvbB02XHLW4HMsafrKAYExjfAl2c4SgRmeB0lPmB0y8bzrNMjj3TJk1EMryW7y8iXzvPrLgpsNGNqR4rKMb4mdPP1nA8dyiDp0v08araD8gHfluNt+mHp/uMU4AHxHn4mzHBxAqKMf4T7zw9cSGwjZ961d2uqt87w+/g6UajKtNVtVBV9wKz+alDRWOCRpTbAYypxwpVdYD3BE8P+P/TPfix8RJ++iMvrpI2lY0b4zrbQzEm8NqJyCnO8BXAd87wFiDVGb6k3DIjRSRORJoDQ/E8mtaYoGIFxZjAWwOMEZHleLoPf9mZ/lfgXyLyLVBabpn5wKfAj3ieax7I55cYUy3W27AxxhifsD0UY4wxPmEFxRhjjE9YQTHGGOMTVlCMMcb4hBUUY4wxPmEFxRhjjE9YQTHGGOMT/w9R1+O30GX5CwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the associated Eq 3.165 Plot\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "pbub = np.linspace(1300,2000)\n",
    "plt.plot(pbub,[Fpbub(pbub, *data) for pbub in pbub])\n",
    "plt.xlabel('Pbub')\n",
    "plt.ylabel('F(Pbub)')\n",
    "plt.title(\"Eq 3.165\")\n",
    "plt.grid();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}