digby1981.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "digby1981.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "authorship_tag": "ABX9TyM390xG92/kYH0y+/hbsek9",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/yohanesnuwara/computational-geophysics/blob/master/digby1981.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_21MieV7pPOf",
        "colab_type": "text"
      },
      "source": [
        "# **Digby (1981) to Model Rock Compressibilities with Respect to Pressure**\n",
        "\n",
        "**Digby (1981) is valid for clastic sedimentary rocks (sandstone, chalk), not for carbonate**\n",
        "\n",
        "Digby (1981) represented a rock as collection of grains with radius `r` and neighboring bonded particles with radius `a`. While **pore pressure** decreases, space `a` decreases to a new radius `b`. Digby (1981) used term **differential pressure** rather than pore pressure. Differential pressure (effective pressure) is confining pressure minus pore pressure: `P_eff = P_conf - Pp`. A new term, coordination number or `Cp` was also introduced. **Murphy (1982)** stated that `Cp` is proportional to `e^(1-porosity)`.  "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Du2GOF94Fq2G",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.optimize import fsolve"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "UhZgwLep1nPW",
        "colab_type": "text"
      },
      "source": [
        "### Step 1. Murphy (1982) to Compute Coordination Number"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gY8kZ4CTts62",
        "colab_type": "text"
      },
      "source": [
        "Murphy (1982) theory of proportionality of coordination number `Cp` and porosity. Supposing the real equation that relates both is `Cp = 11.759 * e^(1 - porosity) - 12.748` [(Li and Ma, 2019)](https://). "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "BwuYwCKEFz67",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 312
        },
        "outputId": "ba3c526c-7b7f-4379-de8a-8ebe5af42114"
      },
      "source": [
        "poro_arr = np.linspace(0, 1, 100)\n",
        "Cp_arr = 11.759 * np.exp(1 - poro_arr) - 12.748\n",
        "plt.plot(poro_arr, Cp_arr)\n",
        "plt.title('Coordination Number to Porosity')\n",
        "plt.xlim(0, 1); plt.ylim(0, 20)\n",
        "plt.xlabel('Porosity'); plt.ylabel('Coordination number')"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "Text(0, 0.5, 'Coordination number')"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 67
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEWCAYAAABIVsEJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0\ndHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd5hU9fXH8fdnG3XpC9KLgDQFAaWr\n2Bt2jR1LwF5iYhKNJsSYn4kmRo2KQUQUFbGLgi0KIl1AuiBI7713OL8/7iVONrvsbJmZ3dnzep77\n7Nx+5u7unLn322RmOOecc/mVkugAnHPOlUyeQJxzzhWIJxDnnHMF4gnEOedcgXgCcc45VyCeQJxz\nzhWIJxAXU5JGS/p5+PoaSZ/H6DwNJO2QlBqL48dT5DVzOQt/100SHUdp5wmkFJB0taQp4T/dakmf\nSOoe7zjM7HUzO7MojiVpiaTTI469zMwqmtnBojh+tnOZpFmSUiKWPSppcFGfKx6yX7sC7G+SdoZ/\nTyslPRnvxB3+rheF8QyW9Gg8z+8CnkCSnKT7gKeA/wNqAQ2A54ELi/g8aUV5vGKoDnBlooPIDwVi\n9T/e1swqAqcBVwN98nuAZLhbLO08gSQxSZWBR4A7zOw9M9tpZvvN7CMzuz/cpoykpyStCqenJJWJ\nOEYfSQslbZI0XFKdiHUm6Q5JC4AF4bIzJM2TtFXSs4Aitr9B0ths+98qaYGkLZKek6Rw3dGSvpK0\nUdIGSa9LqhKuG0KQCD8KvwX/WlKj8Hhp4TZ1wng3hfH3iThvP0lvSXpV0nZJcyR1zONyPg78MadE\nKekUSSuyLfvPt/zwfG9Lei083yxJzSU9IGmdpOWSst+ZHS1psqRtkj6UVC3i2J0ljQ+v2QxJp0Ss\nGy3pz5LGAbuA/3rMk9O1C5dfEF6HLeExWuZxPQAws3nAN0Cb8Dgtw/23hMe7IOLcgyX1lzRS0k6g\np6TK4e9hvaSlkh46nPQkNZX0dfi3tEHSsIhjWbi+L3AN8Ovw/Xwk6X5J72Z7389Iejqa9+Tywcx8\nStIJOBs4AKQdYZtHgIlATSALGA/8KVx3KrABaA+UAf4JjInY14AvgGpAOaAGsB24DEgHfhGe/+fh\n9jcAY7Pt/zFQheBDbT1wdriuKXBGeN4sYAzwVMS+S4DTI+YbhcdLC+fHENxplQXahcc+NVzXD9gD\nnAukAo8BE49wjQxoBkyNeC+PAoPD16cAK7Lt85/4Is53FpAGvAosBn4XXqc+wOKIfUcDKwk+lCsA\n7wKvhevqAhvD2FPCa7QRyIrYdxnQOjxXeg7vJ/u1aw7sDI+VDvwaWAhkHOF6NA1ftwLWADeH+y4E\nHgQyCP5+tgPHhNsOBrYC3cLYy4bX4kMgM/wd/gDcHG4/NLxGh7ftnksMg4FHI9bVDt9PlXA+DVgH\ndEj0/2SyTQkPwKcY/nKDb2Zr8tjmR+DciPmzgCXh65eAxyPWVQT2A43CeSP8UA7nryfig5jg7mMF\nR04gkR8KbwG/zSXOi4DvIuazfwg2Co+XBtQHDgKZEesf46cP/H7AvyPWtQJ2H+EaGUFCOxdYGn44\n5jeBfBGxrhewA0gN5zPDcxz+wBsN/CVbfPsIkt1vgCHZzvUZ0Dti30fy+J1nv3YPA29FzKcQJLBT\njnA9tgGbw7+fR8N9ehAkk5SIbYcC/cLXg4FXI9alhu+rVcSyW4DR4etXgQFAvdx+JxHHfTTb+k+A\nPuHr84G5if5/TMbJH2Elt41AjTzKJ+oQfCgetjRc9j/rzGxHeMy6Edsvz3as5RHbW7b1OVkT8XoX\nQZJCUi1JbyoopN0GvEZwhxONOsAmM9sesWxptrizn7dsHtcJMxtJkBBviTKOSGsjXu8GNthPBf67\nw58VI7aJvG5LCb7d1wAaApeHj4i2SNoCdCf41p3TvtHI/ns+FB6jbq57QHszq2pmR5vZQ+E+dYDl\n4evI2HP7e6kRvq/sf3+Ht/81wZeQyeHjsJvy8Z5eAa4NX18LDMnHvi5KnkCS2wRgL8G399ysIvhQ\nOqxBuOx/1kmqAFQn+HZ6WGR3zqsJvv0f3l6R8/n0f+GxjzWzSgQfAopYf6RupFcB1SRlRixrkC3u\ngvodwSOa8hHLdkbOKygczirkeSKvWwOCO78NBB/AQ8ysSsRUwcz+ErF9Xl1sZ1+f/fd8+PeW3+u1\nCqiv/y64z37dI8+9geB9Zf/7WwlgZmvMrI+Z1SFI2s9LaprDeXN6vx8Ax0lqQ3AH8no+34uLgieQ\nJGZmW4HfA89JukhSeUnpks6R9Hi42VDgIUlZkmqE278Wse5GSe0UFKz/HzDJzJbkcsoRQGtJl4Tf\n5u8Gjipg+JkEj3m2SqoL3J9t/VqyFRAfZmbLCcpyHpNUVtJxBM/oX8tp+/wws9HAbKB3xOIfCO5g\nzpOUDjxEUHZTGNdKaiWpPEE51TvhHctrQC9JZ0lKDd/fKZLq5ePY2a/dW8B5kk4L4/8lwReP8fmM\neRLB3dyvw7+zUwge172Z08bh+3kL+LOkTEkNgfvC94ikyyPe12aCRHEoh0P9z9+Cme0B3gHeACab\n2bJ8vhcXBU8gSc7M/k7wT/kQQUHycuBOgm9oEDy/ngLMBGYB08JlmNm/CZ6Pv0twd3E0R6jKamYb\ngMuBvxA86moGjCtg6H8kKLzfSpCY3su2/jGCxLdF0q9y2P8qgnKRVcD7wB/C91MUHiKoOAD8J1Hf\nDgwk+Pa8k+BRV2EMIXi2v4agAPnu8FzLCapgP8hPv8/7yd//8n9dOzObT3CH90+Cu4JeQC8z25ef\ngMPtewHnhMd5HrjegppaubmL4HotAsYSfOAPCtedAEyStAMYDtxjYduPbF4CWoXv54OI5a8Ax+KP\nr2JGYSGTc84lFUkNgHnAUWa2LdHxJCO/A3HOJZ2wHOY+4E1PHrETswQiqb6kUZLmhjUo7gmXV5P0\nhYLGY19IqprL/r3DbRZI6p3TNs45l11Y2WMbQbuWPyQ4nKQWs0dYkmoDtc1sWlgbZipBbaAbCKpY\n/kXSb4GqZvabbPtWI3gu35Gg4GwqQSOgzTEJ1jnnXL7F7A7EzFab2bTw9Xbge4L63RcSFG4R/syp\niulZBA2vNoVJ4wuCVtXOOeeKibh0gCepEXA8QTW/Wma2Oly1hqCDv+zq8t8NjlaQS6OmsC+cvgAV\nKlTo0KJFi6IJ2jnnSoGpU6duMLMCtVuKeQKRVJGgGui9ZrYtaKMUMDOTVKhnaGY2gKC7Azp27GhT\npkwpzOGcc65UkbQ0761yFtNaWGGjpHeB183scD3+tWH5yOFyknU57LqS/26JW4+iaUXsnHOuiMSy\nFpYIGvh8b2ZPRqwazk+teHsT9MSZ3WfAmZKqhrW0zgyXOeecKyZieQfSDbgOOFXS9HA6l6CV8hkK\nxpA4PZxHUkdJAwHMbBPwJ+DbcHokXOacc66YSKqW6F4G4pxz+SNpqpnlNaBajrwlunPOuQLxBOKc\nc65AkiqBbN9zINEhOOdcqZFUCWTpxp2M/3FDosNwzrlSIakSSEZaCjcPnsK3S7zClnPOxVpSJZAm\nWRWpXaUsN778Ld8t834XnXMulpIqgaSliDd+3plqFTK4ftBkZq/cmuiQnHMuaSVVAgE4qnJZ3ujT\niUpl07lm4CTmrvKxZJxzLhaSLoEA1KtanqF9OlM+I5VrX5rE/DXbEx2Sc84lnaRMIAANqgdJJC1F\nXDNwIgvXeRJxzrmilLQJBKBRjQoM7dsZSVz14iR+XL8j0SE551zSSOoEAnB0VkWG9umEmXHVgIks\n8iTinHNFIukTCEDTmpm80aczBw8ZV704kSUbdiY6JOecK/FKRQIBaF4rk9f7dGL/QU8izjlXFEpN\nAgFocVQlXv95J/bsP+hJxDnnCqlUJRCAlrUr8Uafzp5EnHOukEpdAoH/TiJXDvAk4pxzBVEqEwj8\nlET2HTzElQMmstiTiHPO5UvMEoikQZLWSZodsWxYxPjoSyRNz2XfJZJmhdvFbIzaIIl0CpPIBG8n\n4pxz+RDLO5DBwNmRC8zsZ2bWzszaAe8C7x1h/57htgUaqzdaLY6qxNA+nTlwMGgnsnCdJxHnnItG\nzBKImY0BchyYQ5KAK4ChsTp/fhxzVCZD+3bmkBlXDpjIgrXe7YlzzuUlUWUgPYC1ZrYgl/UGfC5p\nqqS+8Qioea1M3uzbGQmuHDCReWu8F1/nnDuSRCWQqzjy3Ud3M2sPnAPcIemk3DaU1FfSFElT1q9f\nX6igmtbMZFjfzqSnpnDVgInMWeXjiTjnXG7inkAkpQGXAMNy28bMVoY/1wHvAyceYdsBZtbRzDpm\nZWUVOr4mWRUZdktnyqWncvWLk5i5Ykuhj+mcc8koEXcgpwPzzGxFTislVZCUefg1cCYwO6dtY6Vh\n9QoMu6ULlcqlcc2Lk5i61IfHdc657GJZjXcoMAE4RtIKSTeHq64k2+MrSXUkjQxnawFjJc0AJgMj\nzOzTWMWZm/rVyjOsbxdqZJbh+pcmMWnRxniH4JxzxZrMLNExFJmOHTvalClF22xk7bY9XP3iRFZu\n2c3A60+ge7MaRXp855xLJElTC9pcotS2RI9WrUplebNvFxpVr8BNr3zLV/PWJjok55wrFjyBRCEr\nswxD+3TmmFqZ3DJkKp/MWp3okJxzLuE8gUSpaoUMXu/TiePqVeHOod/x/nc51gFwzrlSwxNIPlQq\nm86rN53IiY2qcd9bM3hj0rJEh+SccwnjCSSfKpRJ4+UbT6DnMTV58P1ZDPxmUaJDcs65hPAEUgBl\n01N54doOnHdsbR4d8T1P/fsHkqk2m3PORSMt0QGUVBlpKTx9ZTvKpqfy1L8XsHPvAR48tyVBP5HO\nOZf8PIEUQlpqCk9cdhwVy6Ty4jeL2bH3AI9edCypKZ5EnHPJzxNIIaWkiH4XtCazbDrPjlrIjr0H\nefKKtqSn+tNB51xy8wRSBCTxq7OOoWLZNP7yyTx27NnP89d0oFxGaqJDc865mPGvyUXo1pOP5s8X\nt2H0D+vpPWgy2/bsT3RIzjkXM55Aitg1nRryzJXHM23ZZq5+cSIbd+xNdEjOORcTnkBioFfbOrzY\nuyML1+3g8hcmsHLL7kSH5JxzRc4TSIz0PKYmQ27uxPode7ms/3gWrvNx1p1zycUTSAyd0Kgaw/p2\nYf9B4/IXJjBjuY9u6JxLHp5AYqxVnUq8c2sXKpRJ4+oXJzJ2wYZEh+Scc0XCE0gcNKpRgXdv60r9\nauW5afC3jJjp3cE750o+TyBxUqtSWYbd0oW29Stz59BpDJm4NNEhOedcoXgCiaPK5dIZcnMnTmtR\nk4c/mO2dMDrnSrSYJRBJgyStkzQ7Ylk/SSslTQ+nc3PZ92xJ8yUtlPTbWMWYCId78r2sQz2e+vcC\nHvpgNgcPeRJxzpU8sezKZDDwLPBqtuX/MLO/5baTpFTgOeAMYAXwraThZjY3VoHG2+FOGLMyy9B/\n9I9s3LGPp8KefZ1zrqSI2R2ImY0BNhVg1xOBhWa2yMz2AW8CFxZpcMWAJH5zdgsePr8Vn85ZQ+9B\nk9m627s+cc6VHIkoA7lT0szwEVfVHNbXBZZHzK8Il+VIUl9JUyRNWb9+fVHHGnM3d2/M01e2Y9qy\nzVzxwgTWbN2T6JCccy4q8U4g/YGjgXbAauDvhT2gmQ0ws45m1jErK6uwh0uIC9vVZfCNJ7Jyy24u\neX6ct1p3zpUIcU0gZrbWzA6a2SHgRYLHVdmtBOpHzNcLlyW1bk1rMOyWzuw/ZFzafwLfLinI0z/n\nnIufIyYQSamSRhXVySTVjpi9GJidw2bfAs0kNZaUAVwJDC+qGIqz1nUq895tXaleIYNrBk7i09ne\n4NA5V3wdMYGY2UHgkKTK+T2wpKHABOAYSSsk3Qw8LmmWpJlAT+AX4bZ1JI0Mz3kAuBP4DPgeeMvM\n5uT3/CVV/Wrleee2rrSpU4nbXp/G4HGLEx2Sc87lSHk1ZJP0IXA88AWw8/ByM7s7tqHlX8eOHW3K\nlCmJDqNI7N53kHve/I7P566l70lN+O3ZLUjxsdadc0VM0lQz61iQfaNpB/JeOLk4KpeRSv9rO9Bv\n+BwGjFnEqi27+dvlbb2tiHOu2MgzgZjZK5LKAQ3MbH4cYnKh1BTxyIWtqVe1HI99Mo912/cy4LoO\nVCmfkejQnHMu71pYknoB04FPw/l2kkpFoXZxIIlbTj6aZ646nunLtnBp//Es37Qr0WE551xU1Xj7\nEVS33QJgZtOBJjGMyeXggrZ1GHLziWzYsY+Lnx/ng1M55xIumgSy38y2Zlt2KBbBuCPr1KQ6797W\nlbLpqfxswAQ+n7Mm0SE550qxaBLIHElXA6mSmkn6JzA+xnG5XDStWZH3b+/GMbUyueW1qQwa69V8\nnXOJEU0CuQtoDewFhgLbgHtjGZQ7sqzMMrzZtwtntKzFIx/Ppd/wOd4lvHMu7vJsB/KfDaVKgJlZ\nse2oKZnagUTj4CHjsZHfM3DsYk5vWZOnrzyeCmVi2UO/cy7ZFKYdSDS1sE6QNAuYCcySNENSh4Kc\nzBWt1BTx0Pmt+NOFrflq3jqu+Jf35uuci59oHmG9BNxuZo3MrBFwB/ByTKNy+XJdl0a81PsElmzY\nyYXPjWX2yux1HpxzruhFk0AOmtk3h2fMbCxwIHYhuYLo2aIm79zWlVSJy1+YwBdz1yY6JOdckss1\ngUhqL6k98LWkf0k6RdLJkp4HRsctQhe1lrUr8cEd3WheqyJ9h0xh4DeLiLaMyznn8utIJa7ZB3v6\nQ8Rr/1QqpmpWKsubfbvwy7en8+iI7/lx/Q4eubAN6amJGHzSOZfMck0gZtYznoG4olMuI5Vnr2rP\nkzV+4NlRC1m6cRfPX9Pe+9ByzhWpPOt8SqoCXA80ity+OHbn7n6SkiJ+ddYxNMmqwG/fncXFz49n\nYO+OHJ1VMdGhOeeSRDTPNUYSJI9ZwNSIyZUAl7Svxxt9OrFt934ufm4cYxdsSHRIzrkkEc2AUtPM\nrH2c4imU0taQMD+Wb9rFz1+ZwsL1O+jXqxXXdWmU6JCcc8VATBsSAkMk9ZFUW1K1w1NBTuYSp361\n8rx7e1dOaZ7Fwx/O4aEPZrH/oPeJ6ZwruGgSyD7gCYLxzQ8/vsrza76kQZLWSZodsewJSfMkzZT0\nfli+ktO+S8Kx06dL8luKIlKxTBoDru/ILSc34bWJy7j+pcls3rkv0WE550qoaBLIL4GmYUv0xuEU\nzXggg4Gzsy37AmhjZscBPwAPHGH/nmbWrqC3Vi5nqSnigXNa8vfL2zJ16WYuen4cP6wttt2bOeeK\nsWgSyEIg30PgmdkYYFO2ZZ+b2eFW7BOBevk9risal3aox9C+ndm59yAXPzeOf3vLdedcPkWTQHYC\n08PW6M8cnorg3DcBn+SyzoDPJU2V1PdIB5HUV9IUSVPWr19fBGGVHh0aVuWju7rRJKsifYZM4blR\nC73lunMuatH0/f1BOBUZSb8j6E/r9Vw26W5mKyXVBL6QNC+8o/kfZjYAGABBLayijLM0qF25HG/f\n2oVfvzOTJz6bz/ert/H4ZcdRPsO7hXfOHVmenxJm9kpRnlDSDcD5wGmWy9ddM1sZ/lwn6X2CMdlz\nTCCu8Mqmp/L0le1oWbsSj382j0XrdzLg+g7Uq1o+0aE554qxaMYDWSxpUfapICeTdDbwa+ACM8ux\nXEVSBUmZh18DZwKzc9rWFR1J3HbK0QzqfQLLN+/igmfHMXHRxkSH5ZwrxqIpA+kInBBOPYBngNfy\n2knSUIKqv8dIWiHpZuBZIJPgsdR0SS+E29aRNDLctRYwVtIMYDIwwsw+zef7cgXUs0VNPryjG1XL\np3PtwEm8Mn6Jl4s453IU9ZC2/7VT0HKx2I1K6C3Ri862Pfu5b9h0/v39Oi7rUI9HL2pD2fTURIfl\nnCtihWmJHk1nipHdmKQQ3JF4CWuSq1Q2nQHXdeSpLxfwzJcLWLB2O/2v7UCdKuUSHZpzrpiIJhFE\njgtyAFgCXBGTaFyxkpIi7jujOa3rVOK+YdPp9c+xPHdNezo3qZ7o0JxzxUCBHmEVV/4IK3YWrttO\n3yFTWbpxFw+d15IbujZCUqLDcs4VUqwfYZUBLuV/xwN5pCAndCVT05qZfHBHN+4bNoM/fjSXGcu3\n8Nglx1Euw8tFnCutoqmF9SFwIcHjq50RkytlgnKRDtx3RnM+nLGKS/qPZ9nGfPdy45xLEtGMBzLb\nzNrEKZ5C8UdY8TNq/jrufXM6ZsbTVx5PzxY1Ex2Sc64AYj0eyHhJxxbk4C559TymJh/d2Z16Vctz\n0yvf8o8vfuDQoeQpT3PO5S2aBNIdmCppfjiOxyxJM2MdmCv+GlQvz7u3deWS4+vx9JcLuHHwt2zZ\n5eOLOFdaRFON95yYR+FKrHIZqfzt8uNo37AKfxw+l/P/OZb+13Tg2HqVEx2acy7G8rwDMbOlOU3x\nCM6VDJK4plND3rq1C4cOGZe+MJ43Ji3zLlCcS3LRPMJyLirt6lfh47t70LlJdR58fxa/ensmu/cd\nTHRYzrkY8QTiilS1Chm8fMMJ3Ht6M977bgUXPz+ORet3JDos51wMeAJxRS41Rdx7enMG33gia7ft\n4YJnxzFi5upEh+WcK2LRjAdyiaQFkrZK2iZpu6Rt8QjOlWwnN89ixN09aF6rIne8MY1+w+ew78Ch\nRIflnCsi0dyBPE4wAFRlM6tkZplmVinWgbnkUKdKOd7s24WbujVm8PglXP7CeJZv8tbrziWDaBLI\nWjP7PuaRuKSVkZbC73u14oVr27Now07Oe+Ybvpi7NtFhOecKKZoEMkXSMElXhY+zLpF0Scwjc0nn\n7Da1GXFXDxpUL0+fV6fw5xFz2X/QH2k5V1JFk0AqAbsIxibvFU7nxzIol7wOt16/vktDXvxmMVf8\nawIrNvsjLedKopiOByJpEEGyWXe4Q0ZJ1YBhBN3DLwGuMLPNOezbG3gonH3UzF7J63zemWLJMnLW\nan7zzkxSUsQTlx3Hma2PSnRIzpU6Me1MUVI9Se9LWhdO70qqF+XxBwNnZ1v2W+BLM2sGfBnOZz9n\nNeAPQCfgROAPkqpGeU5XQpx7bG0+vrs79auVo++QqfzxoznsPeAND50rKaJ5hPUyMByoE04fhcvy\nZGZjgE3ZFl8IHL6beAW4KIddzwK+MLNN4d3JF/xvInJJoGH1Crx7W1du6NqIl8ct4bL+E1i60Yeb\nca4kiCaBZJnZy2Z2IJwGA1mFOGctMzvcqmwNUCuHbeoCyyPmV4TL/oekvpKmSJqyfv36QoTlEqVM\nWir9LmjNgOs6sGzTLs57ZiwfTl+Z6LCcc3mIJoFslHStpNRwuhbYWBQnt6AAplCFMGY2wMw6mlnH\nrKzC5DWXaGe2PoqR9/TgmKMyuefN6fz6nRns2ncg0WE553IRTQK5CbiC4G5hNXAZcGMhzrlWUm2A\n8Oe6HLZZCdSPmK8XLnNJrm6Vcgzr25k7ezbl7akr6PXPsXy/2js+cK44irY79wvMLMvMaprZRWa2\nrBDnHA70Dl/3JhhzPbvPgDMlVQ0Lz88Ml7lSIC01hV+ddQyv39yJ7XsOcOFz43hl/BLvHt65YibX\narySfm1mj0v6Jzk8ZjKzu/M8uDQUOAWoAawlqFn1AfAW0ABYSlCNd5OkjsCtZvbzcN+bgAfDQ/3Z\nzPIsuPdqvMln44693P/OTL6at47TW9bk8cvaUq1CRqLDci5pFKYa75ESSC8z+yhsj/E/ommXEW+e\nQJKTmTF4/BIeGzmPKuXT+cfP2tGtaY1Eh+VcUohJOxAz+yh8ucvMXomcCFqmOxcXkrixW2M+uKMb\nlcqlc+1Lk3jsk++9Z1/nEiyaQvQHolzmXEy1qlOJj+7sztUnNuBfXy/i0v7jfbAq5xIoLbcVks4B\nzgXqSnomYlUlwOtWuoQol5HKny8+lpOaZ/Gbd2dy3jNj+UOvVvzshPpISnR4zpUqR7oDWQVMAfYA\nUyOm4QQtxZ1LmLNaH8Wn95xE+4ZV+O17s7j1tals3rkv0WE5V6rk2ZmipHQz2x+neArFC9FLn0OH\njJfGLuaJz+ZTpXw6f7u8LSc19walzkUrpp0pAo0kvSNprqRFh6eCnMy5opaSIvqc1IT37+hK5XLp\nXD9oMn/8aA579nunjM7FWrSdKfYnKPfoCbwKvBbLoJzLr9Z1KvPRXd3/0ynjBc+OZc6qrYkOy7mk\nFk0CKWdmXxI87lpqZv2A82IblnP5VzY96JTxlZtOZPOu/Vz03Dj6j/6Rg4e8BbtzsRBNAtkrKQVY\nIOlOSRcDFWMcl3MFdnLzLD6/9yROb1mLv346j6sGTGT5Jm+65FxRiyaB3AOUB+4GOgDX8VNfVs4V\nS1UrZPD8Ne158oq2fL96G2c/NYa3vl3u/Wk5V4RiOqRtvHktLJeTFZt3cf/bM5mwaCOnt6zFY5cc\nS1ZmmUSH5VyxEOshbZtLelHS55K+OjwV5GTOJUK9quV5/eedePj8VoxZsJ6znhrDJ7NW572jc+6I\nomkHMgN4gaAR4X/qRprZ1NiGln9+B+LysmDtdn7x1nRmr9zGxcfXpV+v1lQun57osJxLmMLcgeTa\nlUmEA2bWvyAHd664aVYrk/dv78Zzoxby7FcLmfDjRv562XGc7I0Pncu3aArRP5J0u6TakqodnmIe\nmXMxkp6awr2nN+f927uRWTaN3oMm88B7M9mx17t4cy4/onmEtTiHxWZmTWITUsH5IyyXX3v2H+Qf\n//6BF8csonblcjxx2XF09bFGXCkS00J0M2ucw1TskodzBVE2PZUHzmnJ27d2JSMthasHTuLhD2az\n0+9GnMvTkbpzP9XMvpJ0SU7rzey92IXlXHx1aFiVkXf34InP5vPy+MWM/mEdj1/ali5HV090aM4V\nW0e6Azk5/Nkrh+n8gp5Q0jGSpkdM2yTdm22bUyRtjdjm9wU9n3PRKpeRyu97teKtW7qQKnHVixP5\n/Yd+N+JcbhLakFBSKrAS6GRmSyOWnwL8yszylai8DMQVld37DvK3z+czaNxi6lYpx18vPc7HYXdJ\nKSbVeCXdd6QdzezJgpwwm9OAHyOTh3PFQbmMVB4+vxXntDmK+9+ZyTUDJ3HViQ144NwWVCrr7Uac\ngyM/wsoMp47AbUDdcLoVaMoG/+EAABW7SURBVF9E578SGJrLui6SZkj6RFLr3A4gqa+kKZKmrF+/\nvojCci7QsVE1PrmnB31PasKwb5dx1j/GMGr+ukSH5VyxEE013jHAeWa2PZzPBEaY2UmFOrGUQTBs\nbmszW5ttXSXgkJntkHQu8LSZNcvrmP4Iy8XS9OVb+PU7M/hh7Q4uOb4uD5/fiqoVMhIdlnOFEusR\nCWsBkYNN7wuXFdY5wLTsyQPAzLaZ2Y7w9UggXZI/gHYJ1a5+FT66qzt3ndqU4TNWccY/vmbkrNXe\nw68rtaJJIK8CkyX1k9QPmAQMLoJzX0Uuj68kHSVJ4esTwzg3FsE5nSuUMmmp/PLMYxh+Z3dqVy7H\n7a9P49bXprJ2255Eh+Zc3EVVC0tSe6BHODvGzL4r1EmlCsAyoImZbQ2X3QpgZi9IupOg3OUAsBu4\nz8zG53Vcf4Tl4unAwUO8NHYxT37xAxlpKTx4bkuuPKE+4Xcf50qEwjzCOmICCavZzjGzFgUNLp48\ngbhEWLxhJw+8N5OJizbRuUk1HrvkOBrXqJDosJyLSszKQMzsIDBfUoMCReZcKdC4RgWG9unMXy45\nljmrtnHWU2N4btRC9h88lOjQnIupaLpzrwrMkTQZ2Hl4oZldELOonCthJHHliQ04tUVN+n00hyc+\nm89HM1bxl0uPo139KokOz7mYiKYa78k5LTezr2MSUSH4IyxXXHwxdy0PfzCbtdv30LtLI355ZnMy\nvQGiK4Zi3Rvv18A8fmpY+H1xTB7OFSdntKrFF/edRO8ujXhlwhLOeHIMn81Zk+iwnCtS0YyJfgUw\nGbgcuAKYJOmyWAfmXEmXWTadfhe05r3bulKlfDq3DJlKn1ensGrL7kSH5lyRiHZM9DPMbF04nwX8\n28zaxiG+fPFHWK642h9W+X3q3z+QInHfGc25oWsj0lKjaYrlXOzEuiV6yuHkEdoY5X7OuVB6agq3\nnnw0X/ziZDo1rsajI77nwufGMX35lkSH5lyBRZMIPpX0maQbJN0AjABGxjYs55JT/WrlGXTDCTx3\ndXs27NjLxc+P4+EPZrN19/5Eh+ZcvkXbEv0SoHs4+42ZvR/TqArIH2G5kmT7nv08+cUPvDJ+CdUq\nlOHh81tyQds63pLdxVWsH2EBjANGAV+Fr51zhZRZNp0/9GrN8Du7U7dKWe55czpXvziJhet2JDo0\n56KSn1pYl+G1sJwrcm3qVua927vx6EVtmLNqK+c8PYbHP53Hrn0+lK4r3rwWlnPFyIYde3ls5Dze\nnbaCulXK8fterTizVS1/rOVixmthOZckalQsw9+vaMtbt3Qhs2watwyZyo2Dv2XJhp157+xcnBW0\nFtYnsQ3LudLtxMbV+Oiu7jx0XkumLNnMmf8Yw98/n8/ufQcTHZpz/+G1sJwr5tZt28P/jfyeD6av\nom6Vcjx8fivOau2PtVzRiMl4IJKaArXMbFy25d2B1Wb2Y0FOGEueQFwym7RoI7//cA7z126nR7Ma\n9LugNUdnVUx0WK6Ei1UZyFPAthyWbw3XOefiqFOT6oy4uzv9erVi+vItnP3UGB4b+T3b93gjRJcY\nR0ogtcxsVvaF4bJGMYvIOZertNQUbujWmFG/OoWLj6/Lv8Ys4tS/f827U1dw6FDej6OdK0pHSiBH\nGgWnXGFPLGmJpFmSpkv6n+dOCjwjaaGkmeG47M45gtpaj1/Wlg/u6EadKuX45dszuPSF8czwvrVc\nHB0pgUyR1Cf7Qkk/B6YW0fl7mlm7XJ6/nQM0C6e+QP8iOqdzSaNd/Sq8f1tXnrjsOJZv2s2Fz43j\n/rdnsG77nkSH5kqBIw1pey/wvqRr+ClhdAQygItjHRhwIfCqBaX8EyVVkVTbzFbH4dzOlRgpKeLy\njvU5u81RPPvVQgaNW8wns9dw56lNubFbI8qkpSY6RJekcr0DMbO1ZtYV+COwJJz+aGZdzKwohlYz\n4HNJUyX1zWF9XWB5xPyKcNl/kdRX0hRJU9avX18EYTlXMmWWTeeBc1vyedhl/F8+mceZ/xjD53PW\nEE11fefyK6p2IDE5sVTXzFZKqgl8AdxlZmMi1n8M/MXMxobzXwK/MbNc6+l6NV7nfjLmh/X86eO5\nLFi3g25Nq/PQea1oWbtSosNyxUw8euMtcma2Mvy5DngfODHbJiuB+hHz9cJlzrkonNQ8i5H39OCP\nF7RmzqptnPfMNzzw3iw27Nib6NBckkhIApFUQVLm4dfAmcDsbJsNB64Pa2N1BrZ6+Ydz+ZOemkLv\nro0Y/atT6N21EW9PWc4pT4ym/+gf2bPfu0VxhZOoO5BawNiwp9/JwAgz+1TSrZJuDbcZCSwCFgIv\nArcnJlTnSr4q5TP4Q6/WfPaLk+jUuBp//XQepz/5NR/PXOXlI67AElYGEgteBuJcdMYt3MCfPp7L\nvDXbad+gCr87rxUdGlZNdFguAUpkGYhzLnG6Na3BiLt78NdLj2X55t1c2n88d7wxjeWbdiU6NFeC\n+B2Ic6Xczr0H+NeYRQwY8yOHDkHvrg25s2czKpdPT3RoLg78DsQ5V2AVyqRx3xnNGf2rnlx0fB0G\njl3MSU+MYuA3i9h7wAvaXe48gTjnADiqclkev6wtI+/uQdv6VXh0xPec9vev+XD6Su+o0eXIE4hz\n7r+0rF2JV286kVdvOpHMsunc8+Z0LnxuHOMXbkh0aK6Y8QTinMvRSc2zGHFXd568oi2bdu7j6oGT\n6D1oMnNX5TRMkCuNPIE453KVkiIuaV+PL395Mg+e24Lpy7dw3j+/4b5h073GlvNaWM656G3dtZ/+\nX//Iy+MWYwbXdG7AnT2bUr1imUSH5gooJmOil0SeQJyLj9Vbd/PUFwt4e+pyymek0adHE27u0ZiK\nZY40QoQrjjyBhDyBOBdfC9dt54nP5vPZnLVUr5DBnac25epODXwMkhLEE0jIE4hzifHdss08/ul8\nJizaSN0q5bj39GZcfHxd0lK9mLW484aEzrmEOr5BVd7o04khN59I9YoZ3P/OTM56agwjZ632zhqT\nmCcQ51yRkESPZll8eEc3Xri2PZK4/fVp9Hp2LKPmr/NEkoQ8gTjnipQkzm5Tm8/uPYm/Xd6WLbv2\nc+PL33LFvyYwcdHGRIfnipCXgTjnYmrfgUMMm7Kcf365gHXb99K9aQ3uO7M57Rt49/HFgReihzyB\nOFd87dl/kNcmLqX/6B/ZuHMfp7aoyX1nNKdN3cqJDq1U8wQS8gTiXPG3c+8BBo9fwoAxi9i6ez9n\ntqrFL85oTsvalRIdWqnkCSTkCcS5kmPbnv28PHYJA79ZxPa9Bzj32KO457TmHHNUZqJDK1VKVDVe\nSfUljZI0V9IcSffksM0pkrZKmh5Ov493nM652KpUNp17Tm/G2N+cyl2nNmXMDxs4++kx3PHGNBas\n3Z7o8FwU4n4HIqk2UNvMpknKBKYCF5nZ3IhtTgF+ZWbn5+fYfgfiXMm1eec+Bo5dxOBxS9i1/yDn\nHVubu09rRvNafkcSSyXqDsTMVpvZtPD1duB7oG6843DOFS9VK2Rw/1ktGPubU7n9lKMZNW8dZz0V\n3JHMX+N3JMVRQtuBSGoEHA9MymF1F0kzJH0iqXVcA3POJUz2RDI6TCS3vz6V71f7WCTFScIK0SVV\nBL4G/mxm72VbVwk4ZGY7JJ0LPG1mzXI5Tl+gL0CDBg06LF26NMaRO+fiafPOfQwat5jB45awfe8B\nzmxVi7tPa+bVf4tIiauFJSkd+Bj4zMyejGL7JUBHMzvimJpeBuJc8tq6az+Dxi3m5XGL2bbnAD2P\nyeKu05p5g8RCKlEJRJKAV4BNZnZvLtscBaw1M5N0IvAO0NDyCNYTiHPJb9ue/bw6fgkDxy5my679\ndG9agztPbUqnxtUIPl5cfpS0BNId+AaYBRwKFz8INAAwsxck3QncBhwAdgP3mdn4vI7tCcS50mPn\n3gO8PmkpA8YsZsOOvZzQqCp39GzKyc2zPJHkQ4lKILHkCcS50mfP/oMM+3Y5//r6R1Zt3UObupW4\n45SmnNX6KFJSPJHkxRNIyBOIc6XXvgOH+OC7lfT/+kcWb9hJ05oVufXko7mwXR3SfWCrXHkCCXkC\ncc4dPGSMmLWa50ctZN6a7dStUo6+JzXhio71KZfhQ+1m5wkk5AnEOXeYmTF6/nqeG7WQKUs3U71C\nBjd2a8R1XRpRuVx6osMrNjyBhDyBOOdyMnnxJvqPXsio+eupWCaNqzs14ObujalVqWyiQ0s4TyAh\nTyDOuSOZu2obL3z9Ix/PXEVqirj4+Lr0PelomtasmOjQEsYTSMgTiHMuGss37eLFbxbx1pTl7D1w\niNNb1uLWk5vQoWG1RIcWd55AQp5AnHP5sXHHXl6ZsJRXJyxhy679dGhYlb4nNeGMlrVKTRVgTyAh\nTyDOuYLYte8Ab09ZwYvfLGLF5t00rlGBn/dozKXt61E2PblrbnkCCXkCcc4VxoGDh/hk9hoGjFnE\nrJVbqV4hg+u6NOS6zg2pXrFMosOLCU8gIU8gzrmiYGZMWLSRgd8s5qt56yiTlsIl7evx8x6NOTor\nuQrcC5NA0oo6GOecK+kk0fXoGnQ9ugYL121n4DeLeXfaCoZOXsZpLWpyc4/GdGlSvdT3ueV3IM45\nF4X12/fy2sSlDJm4lE0799GqdiVu7t6YXm3rkJFWcrtK8UdYIU8gzrlY27P/IO9/t5JBYxezYN0O\nsjLLcH3nhlzdqUGJLCfxBBLyBOKcixczY8yCDbw0djFjflhPmbQULj6+Ljd2a8wxR2UmOryoeRmI\nc87FmSRObp7Fyc2zWLB2Oy+PX8J701bw5rfL6da0Ojd2bcypLWomdXsSvwNxzrkisnnnPoZ+u4wh\nE5ayeuseGlYvT+8ujbisYz0qlS2eHTj6I6yQJxDnXHGw/+AhPpuzhpfHLWHq0s1UyEjl0g71uL5L\no2LX75YnkJAnEOdccTNrxVYGj1/CRzNWse/gIXo0q0HvLo3o2aImqcXg8ZYnkJAnEOdccbVhx17e\nnLyM1yYuY822PdSvVo7rOjfkio71qVI+I2FxlbgEIuls4GkgFRhoZn/Jtr4M8CrQAdgI/MzMluR1\nXE8gzrnibv/BQ3w+Zy2vjF/C5CWbKJOWwoXt6nB9l0a0qVs57vGUqAQiKRX4ATgDWAF8C1xlZnMj\ntrkdOM7MbpV0JXCxmf0sr2N7AnHOlSTz1mzj1QlLeX/aSnbvP8jxDapwfZeGnNOmdtw6cSxpCaQL\n0M/MzgrnHwAws8citvks3GaCpDRgDZBleQTrCcQ5VxJt3b2f96atYMjEpSxav5NqFTK4vGM9rjmx\nIQ2ql4/puUtaO5C6wPKI+RVAp9y2MbMDkrYC1YEN2Q8mqS/QN5zdK2l2kUdcMtUgh+tVCvl1+Ilf\ni58U62uxFPgOeDA+pzumoDuW+IaEZjYAGAAgaUpBM2my8WsR8OvwE78WP/Fr8RNJBX5sk4gewFYC\n9SPm64XLctwmfIRVmaAw3TnnXDGRiATyLdBMUmNJGcCVwPBs2wwHeoevLwO+yqv8wznnXHzF/RFW\nWKZxJ/AZQTXeQWY2R9IjwBQzGw68BAyRtBDYRJBkojEgJkGXTH4tAn4dfuLX4id+LX5S4GuRVA0J\nnXPOxU/JHQXFOedcQnkCcc45VyAlLoFIOlvSfEkLJf02h/VlJA0L10+S1Cj+UcZHFNfiPklzJc2U\n9KWkhomIMx7yuhYR210qySQlbRXOaK6FpCvCv405kt6Id4zxEsX/SANJoyR9F/6fnJuIOONB0iBJ\n63JrK6fAM+G1mimpfZ4HNbMSMxEUuv8INAEygBlAq2zb3A68EL6+EhiW6LgTeC16AuXD17eV5msR\nbpcJjAEmAh0THXcC/y6aEbRTqxrO10x03Am8FgOA28LXrYAliY47htfjJKA9MDuX9ecCnwACOgOT\n8jpmSbsDORFYaGaLzGwf8CZwYbZtLgReCV+/A5wmKfF9Jhe9PK+FmY0ys13h7ESCNjfJKJq/C4A/\nAX8F9sQzuDiL5lr0AZ4zs80AZrYuzjHGSzTXwoBK4evKwKo4xhdXZjaGoFZrbi4EXrXARKCKpNpH\nOmZJSyA5dYNSN7dtzOwAcLgblGQTzbWIdDPBt4tklOe1CG/H65vZiHgGlgDR/F00B5pLGidpYtg7\ndjKK5lr0A66VtAIYCdwVn9CKpfx+ppT8rkxc3iRdC3QETk50LIkgKQV4ErghwaEUF2kEj7FOIbgr\nHSPpWDPbktCoEuMqYLCZ/T3s6HWIpDZmdijRgZUEJe0OxLtB+Uk01wJJpwO/Ay4ws71xii3e8roW\nmUAbYLSkJQTPd4cnaUF6NH8XK4DhZrbfzBYTDK/QLE7xxVM01+Jm4C0AM5sAlCXoaLE0iuozJVJJ\nSyDeDcpP8rwWko4H/kWQPJL1OTfkcS3MbKuZ1TCzRmbWiKA86AIzS8a+/6P5H/mA4O4DSTUIHmkt\nimeQcRLNtVgGnAYgqSVBAlkf1yiLj+HA9WFtrM7AVjNbfaQdStQjLIttNyglSpTX4gmgIvB2WI9g\nmZldkLCgYyTKa1EqRHktPgPOlDQXOAjcb2ZJd5ce5bX4JfCipF8QFKjfkKRfOJE0lOCLQ42wzOcP\nQDqAmb1AUAZ0LrAQ2AXcmOcxk/RaOeeci7GS9gjLOedcMeEJxDnnXIF4AnHOOVcgnkCcc84ViCcQ\n55xzBeIJxLlcSDooabqk2ZLellQ+hud6JGz0iaR7Y3ku54qKV+N1LheSdphZxfD168BUM3syiv1E\n8L9VoO4wwtbyHc1sQ0H2dy5e/A7Eueh8AzSF/4yzMjuc7g2XNQrHnXgVmA3Ul3SVpFnhdn8Nt0uV\nNDhcNitswEa47DJJdwN1gFHhOBU3SXrqcBCS+kj6R5zfu3M5KlEt0Z1LhLBPtXOATyV1IGih24lg\n3IRJkr4GNhP0J9XbzCZKqkPQdXyHcN3nki4i6O20rpm1CY9dJfJcZvaMpPuAnma2QVJF4HeS7jez\n/eG5b4nD23YuT34H4lzuykmaDkwh6DPpJaA78L6Z7TSzHcB7QI9w+6XhOAoAJwCjzWx9OKzA6wQD\n+iwCmkj6Z9iN+rYjBRCe4yvgfEktgHQzm1W0b9O5gvE7EOdyt9vM2kUuyGNssp15HdDMNktqC5wF\n3ApcAdyUx24DgQeBecDLeZ3DuXjxOxDn8ucb4CJJ5SVVAC4Ol2U3GThZUg1JqQTjTnwd9n6bYmbv\nAg8RDDGa3XaCLugBMLNJBN1sXw0MLdJ341wh+B2Ic/lgZtMkDSZIEAADzew7SY2ybbda0m+BUQRl\nJSPM7MPw7uPlcJArgAdyOM0AgvKWVWbWM1z2FtDu8DC0zhUHXo3XuRJA0sfAP8zsy0TH4txh/gjL\nuWJMUhVJPxCUx3jycMWK34E455wrEL8Dcc45VyCeQJxzzhWIJxDnnHMF4gnEOedcgXgCcc45VyD/\nD84wuvdY+eu0AAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vIitFe0L1w2I",
        "colab_type": "text"
      },
      "source": [
        "### Step 2. Calculate Rock Properties (Mineral Bulk and Shear Incompressibility)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "DHK0qubIvuL7",
        "colab_type": "text"
      },
      "source": [
        "Supposing a sandstone has **50% porosity**, composed of **50% quartz** and **50% clay minerals**, and **Poisson's ratio of 0.4** (Typically sandstone has `poisson` between 0.1 and 0.4). Mineral bulk moduli (`Km`) and shear moduli (`Gm`) can be computed using **Voigt-Reuss-Hill method**. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "cFBjmP8gwJcK",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 67
        },
        "outputId": "9a7d7d14-a8aa-49d4-81fe-4137d5076b09"
      },
      "source": [
        "# Voigt Reuss Hill to compute Km and Gm\n",
        "\n",
        "rhocalc = 2.71; Kcalc = 76.8; Gcalc = 32\n",
        "rhoclay = 2.58; Kclay = 20.9; Gclay = 6.9\n",
        "rhodolo = 2.87; Kdolo = 94.9; Gdolo = 45\n",
        "rhoqtz = 2.65; Kqtz = 36.6; Gqtz = 45\n",
        "# mineral composition\n",
        "calc = 0; clay = 0.5; dolo = 0; qtz = 0.5\n",
        "\n",
        "mincomposition = np.array([calc, clay, dolo, qtz])\n",
        "rhomineral = np.array([rhocalc, rhoclay, rhodolo, rhoqtz])\n",
        "Kmineral = np.array([Kcalc, Kclay, Kdolo, Kqtz])\n",
        "Kmineral_inv = 1 / Kmineral\n",
        "Gmineral = np.array([Gcalc, Gclay, Gdolo, Kqtz])\n",
        "Gmineral_inv = 1 / Gmineral\n",
        "\n",
        "Kv = sum(mincomposition * Kmineral)\n",
        "Kr_inv = sum(mincomposition * Kmineral_inv)\n",
        "Gv = sum(mincomposition * Gmineral)\n",
        "Gr_inv = sum(mincomposition * Gmineral_inv)\n",
        "rhom = sum(mincomposition * rhomineral)\n",
        "Km = (Kv+(1/Kr_inv)) / 2\n",
        "Gm = (Gv+(1/Gr_inv)) / 2\n",
        "\n",
        "print(\"Mineral bulk moduli:\", Km, \"GPa\")\n",
        "print(\"Mineral shear moduli:\", Gm, \"GPa\")\n",
        "print(\"Mineral density:\", rhom, \"g/cc\")"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Mineral bulk moduli: 27.678304347826085 GPa\n",
            "Mineral shear moduli: 16.68051724137931 GPa\n",
            "Mineral density: 2.615 g/cc\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "nVvYAReUIbdx",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        },
        "outputId": "9a885846-2d0e-4b56-e896-afe483a0f1a5"
      },
      "source": [
        "poro = 0.5\n",
        "Cp = 11.759 * np.exp(1 - poro) - 12.748\n",
        "print(\"Rock coordination number:\", Cp)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Rock coordination number: 6.639313422162809\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "U732fLsRxI2Y",
        "colab_type": "text"
      },
      "source": [
        "The sandstone has mineral bulk moduli (`Km`) 27.68 GPa, shear moduli (`Gm`) 16.68 GPa, and coordination number (`Cp`) 6.64. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "I-DWfDPd15R8",
        "colab_type": "text"
      },
      "source": [
        "### Step 3. Calculate Differential Pressure\n",
        "\n",
        "Assuming the sandstone is confined to a pressure of `P_conf = 70 MPa` and the pore pressure is `Pp = 60 MPa`.  "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "WVXJdUrfyWXI",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        },
        "outputId": "8b544895-7eb7-4a81-bf6a-efb85623190b"
      },
      "source": [
        "P_conf = 70 \n",
        "Pp = 60\n",
        "pressure = P_conf - Pp\n",
        "print(\"Differential pressure:\", pressure, \"MPa\")"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Differential pressure: 10 MPa\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "C2s8-GOu2DkI",
        "colab_type": "text"
      },
      "source": [
        "### Step 4. Use Digby (1981) to Compute Dry Rock Incompressibility"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "04N58Agpyl8h",
        "colab_type": "text"
      },
      "source": [
        "Since radius `a` and `r` not known, we assume the rock has ratio of a and r (`ratio_ar`) is 1. We then compute dry rock incompressibility (`K_dry`) using **Digby (1981) method**. \n",
        "\n",
        "First, solve `d`. Second, calculate ratio of deformed radius `b` to `r` (`ratio_br`). Third, compute `K_dry`."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "QcTIptKkGetp",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        },
        "outputId": "24ac17db-18cf-4e05-dd37-c63b45a0084f"
      },
      "source": [
        "\"Digby (1981) Method\"\n",
        "\n",
        "# ratio_ar = np.linspace(0, 1, 100)\n",
        "ratio_ar = 1\n",
        "\n",
        "poisson = 0.4 # typical sandstone, high poisson ratio\n",
        "Gm = 10 # shear moduli, in GPa\n",
        "pressure = 10 # pressure, in MPa\n",
        "Cp = 10\n",
        "\n",
        "def f(d):\n",
        "  f1 = (d**3) + (1.5 * (ratio_ar**2) * d) - ((3 * np.pi * (1 - poisson) * pressure) / (2 * Cp * (1 - poro) * Gm))\n",
        "  return(f1)\n",
        "\n",
        "d = fsolve(f, 1)\n",
        "\n",
        "ratio_br = np.sqrt((d**2) + (ratio_ar**2))\n",
        "Kdry = ((Cp * (1 - poro) * Gm) / (3 * np.pi * (1 - poisson))) * ratio_br\n",
        "print(\"Dry bulk modulus:\", Kdry, \"GPa\")"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Dry bulk modulus: [9.36413669] GPa\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Oz4QUVZs1lK9",
        "colab_type": "text"
      },
      "source": [
        "### Step 5. Zimmermann (1986) to Compute Pore Compressibility\n",
        "\n",
        "Use [(Zimmermann 1986)](https://() to compute pore incompressibility (`K_pore`) from the resulted `K_dry`. Then, invert `K_pore` to arrive at value of **pore compressibility** or the **formation compressibility** (`cf`). \n",
        "\n",
        "`cf` later is important in **reservoir engineering** since it determine the term **total compressibility as sum of formation compressibility and fluid compressibility** or `ct = cf + (sw * cw) + (so * co) + (sg * cg)`, although often negligible. \n",
        "\n",
        "\n",
        "\n",
        "```\n",
        "sw, so, sg = water, oil, gas saturation\n",
        "sw + so + sg = 1\n",
        "cw, co, cg = water, oil, gas isothermal compressibility\n",
        "```\n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "vRo77pIZgYqJ",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        },
        "outputId": "b7eccaf5-f1de-4079-9361-fcba634769a3"
      },
      "source": [
        "K_pore_inverse = (1 / (poro * Kdry)) - ((1 - poro) / (poro * Km))\n",
        "cf = (1 / K_pore_inverse) / 145038 # from GPa^-1 to psi^-1\n",
        "print(\"Pore compressibility:\", cf, \"psi^-1\")"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Pore compressibility: [3.88542625e-05] psi^-1\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8nlvRM9U6R3O",
        "colab_type": "text"
      },
      "source": [
        "### Plot Change of Pore Compressibilities with Declining Pore Pressure"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "11rRPrRi6gFW",
        "colab_type": "text"
      },
      "source": [
        "During production, assuming **pore pressure declines from 60 to 5 MPa**, at **confining pressure 70 MPa**. Sandstone with the same properties, plot the pore compressibility change with respect to declining pore pressure."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "XkdQBdrJ64UR",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 351
        },
        "outputId": "481f615e-7914-4ad1-99e7-a2d314a980c9"
      },
      "source": [
        "P_conf - 70\n",
        "Pp = np.arange(5, 61, 1)\n",
        "pressure = P_conf - Pp\n",
        "\n",
        "cf_arr = []\n",
        "\n",
        "for i in range(len(pressure)):\n",
        "\n",
        "  def f(d):\n",
        "    f1 = (d**3) + (1.5 * (ratio_ar**2) * d) - ((3 * np.pi * (1 - poisson) * pressure[i]) / (2 * Cp * (1 - poro) * Gm))\n",
        "    return(f1)\n",
        "  d = fsolve(f, 1)\n",
        "  \n",
        "  ratio_br = np.sqrt((d**2) + (ratio_ar**2))\n",
        "  Kdry = ((Cp * (1 - poro) * Gm) / (3 * np.pi * (1 - poisson))) * ratio_br\n",
        "\n",
        "  K_pore_inverse = (1 / (poro * Kdry)) - ((1 - poro) / (poro * Km))\n",
        "  cf = (1 / K_pore_inverse) / 145038 # from GPa^-1 to psi^-1\n",
        "  cf_arr.append(float(cf))\n",
        "\n",
        "# plot\n",
        "plt.figure(figsize=(15, 5))\n",
        "\n",
        "plt.subplot(1, 2, 1)\n",
        "plt.plot(pressure, cf_arr, '.-', color='blue')\n",
        "plt.title('Pore Compressibility to Differential Pressure')\n",
        "# plt.xlim(0, 1); plt.ylim(0, 20)\n",
        "plt.xlabel('Differential pressure (MPa)'); plt.ylabel('Pore compressibility (psi^-1)')\n",
        "\n",
        "plt.subplot(1, 2, 2)\n",
        "plt.plot(Pp, cf_arr, '.-', color='red')\n",
        "plt.title('Pore Compressibility to Pore Pressure')\n",
        "# plt.xlim(0, 1); plt.ylim(0, 20)\n",
        "plt.xlabel('Pore pressure (MPa)'); plt.ylabel('Pore compressibility (psi^-1)')\n",
        "\n",
        "plt.tight_layout(pad=2.0) # set distance between two plots"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAABB4AAAFOCAYAAAAhEN5jAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0\ndHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd5xU5fXH8c9hKUssqEgKxoKiRiyx\nbFSsgIZiwY5YgijWWJLYECuKgqjB3nuJBbGhMSLSRF3LYokiGpGYoCaKKIpKWdjz++Pc+e0w2Z2d\nXXZ2tnzfr9e8mHvnzjPPlJ05nPs85zF3R0REREREREQkH1oVugMiIiIiIiIi0nwp8SAiIiIiIiIi\neaPEg4iIiIiIiIjkjRIPIiIiIiIiIpI3SjyIiIiIiIiISN4o8SAiIiIiIiIieaPEg0gLZGbrmdn3\nZlaUbE81s2Nre6yZHWFmzzdcz1dOFc/lZ2b2opktNLM/W7jbzL4xs9cbW39rOLaHmX3aEP0SERGp\nTy01LhFpSZR4kEbFzD4xs0XJD8oXZnaPma3aQI+9upldY2b/Th7/42R77YZ4/Ibk7v9291XdffnK\nHOvuf3H33qltM3Mz61qXPq3sf5zNbLCZLU/eu+/N7J9JEmGTtP5mPpfjga+A1d39DGAX4LfAL919\n+7r2ZSWewydmtmeW/q5M225mPySvzWdmNiaXhIaISEumuKRhtIC45Dsze9vM9qlre7V87Klmtjh5\n7K/M7HEz+0VDPLZIdZR4kMZoX3dfFdgWKAHOr20DZta6lse3BSYBmwN9gdWB7sB8oMH/A5quts+l\nhStNPjsdgD2BRcAMM9uimuPXB953d0/b/sTdf6jtAzeR9+nXyeuzB3A4cFzmAY3heTSGPoiIpFFc\nkkbf0bWSikvWAO4ExprZmrVpYCVe71OSx94kefyrq2m/oCch9HlqOZR4kEbL3T8D/gZsAWBmnc1s\nvJl9bWazzez//9NkZsPNbJyZPWBm3wGDzayVmZ2TnCGYb2ZjzWytah5uELAecIC7v+/uFe7+pbuP\ncPdnk8fYLMkgLzCzmWbWP+3x7zGzm8zsb0l2+WUz+3lyZuIbM/vAzLZJO/4TMxtmZu8nt99tZsXJ\nbT3M7FMzG2pm/wXuTvbvk2TLF5jZK2a2VVp7Q5Oz2AvN7EMz2yPZv72ZlSWZ9i/MbEyyf4PkLED6\nl/1GZvZ6cuxTqdeqmmNTjzvYzF5Krr+Y7H4neQ0ONbP3zGzftOPbJJn3bTLaWSV5rztb5YiFzmbW\nLnkNP08u15hZu2rew//n7svd/WN3/z0wDRie+VzM7B7gKODs5PFOAO4AuifbF+fwun+SvPZ/B35I\n2u1sZo+Z2TyLURenpR0/PPkc3pe8VzPNrCS57X7iM/h08vhnZ772Zna0mc1K7jsn6XOtufsHwHQq\n/7Zq+zyq+1wVW/wNzk9erzfM7Gdpj7FnWhvDzeyBjPdliJn9G5ic7N8xec0XmNk7ZtajLs9XRKQ+\nKC5RXGJ1j0sqgLuA9sBGyWMcl3xuvk4+R53THt/N7GQz+wj4KNn3KzObmBz/oZkNqOlxk8f+GniM\nys/tPWZ2s5k9a2Y/AD2T53WVxeiaL8zsFjNrnxy/tpk9k7zPX5vZdDNrldxW3ft8j5ldmvZ8Vhg9\nYrWMO6SZcHdddGk0F+ATYM/k+rrATGBEsv0icBNQDGwNzAN6JbcNB8qB/YmEWnvgD8CrwC+BdsCt\nwEPVPO7DwL1Z+tUGmA2cC7QFegELgU2T2+8hhuxvl/RvMvBPInAoAi4FpmQ8z/eS57gW8DJwaXJb\nD2AZMDrpd3tgG+BLYIekvaOSNtoBmwJzgc7J/TcANkqulwK/S66vCuyYdowDrZPtqcBnxI/SKsQP\n1ANZjj02uT4YeCnteTnQNW37bOCRtO39gHereY17AJ9m7LskeQ9/CnQCXkl9Hqq4/wp9Sdt/DPBF\nNc/lntTrXs3zqfZ1T3sf307ex/bEZ28GcCHxOdkQmAP0SfucLgb2StobBbxa1ee/mv7uTQQsBuwO\n/AhsW93rl/E6/P97A3QD/gsMqePzqO5zdQLwNPCT5PltR0xjqeq5Ded/P2P3EZ+/9sA6xJm9vZL+\n/DbZ7lTo7ylddNGl5VxQXJL6fVFc4nWPS4DWyfu/kBiV2St5f7ZNXrPrgRcz+j0xeS/aJ6/BXODo\npK1tkvt3q+ax01+TtZP3//60z8a3wM7EZ7OYGA0xPnm81Yjf8lHJ8aOAW4jPXBtgVyIOyfY+38OK\n8dUKryW1jDt0aR4XjXiQxuhJM1sAvEScrR5pZusSX5BD3X2xu79NnJ0elHa/Und/0uOswCLgROA8\nd//U3ZcQQcDBVWXIgY7Af7L0aUfiB/Jyd1/q7pOBZ4DD0o55wt1nuPti4Algsbvf5zEH8RHiRyLd\nDe4+1yMTfVlGWxXARe6+JHkuxwO3uvtrHmfz7wWWJP1aTvxodTOzNu7+ibt/nLRTDnQ1s7Xd/Xt3\nfzXLc7zf3d/zmGZwATDAVn743QPAXma2erL9O+D+Wtz/COASj7M884CLkzZq43Pih7Qusr3uKdcl\n7+Mi4DfEf4wvST4nc4DbgYFpx7/k7s8mn4v7gV/n2hl3/6vHSA5392nA80QAkKs3zewbIqC4g+Ss\nVR2eR3Wfq3Lib6lr8nrNcPfvatG/4e7+Q9KHI4Fnk9eqwt0nAmVEIkJEpCEpLlFcklLbuGTH5LPz\nX+L1PMDdv03aucvd30w+C8OIEZcbpN13lLt/nbze+xBTQe9292Xu/haRjDkky2Nflzz2O8Rn6fS0\n255y95c9RmIsId7PPyWPtxAYyYq/+b8A1nf3cnef7u5O9vc5F7WNn6SJU+JBGqP93X0Nd1/f3X+f\nfCF1BlJfhin/Is6KpszNaGd94IlkaNgCYBbxJfmzKh5zPvGlWp3OwNzkC7q6x/8i7fqiKrYzi1Gl\n9/dfyWOkzEsChfTnckbquSTPZ10iyzwb+CMRwHxpZg+nDdcbQszt+8Bi2Hu2okaZ/WlDZMnrzN0/\nJ86aHGRmawD9gL/UoonOSV/S+9W5mmOrsw7wdS3vk1Lt6552zNyM4ztnHH8uK37m/pt2/UeguJqg\n83+YWT8zezUZ6riA+E94bd6jbd19TXffyN3Pz/g81+Z5VPe5uh+YADycDEG9wsza1KJ/mX04JKMP\nu5D971REJB8UlyguSaltXPJq8tlZ2913dPcXqmrH3b8n3vPqPj/rAztkvN5HAD/P8tinJY+9jrsf\nkSRKqmq7EzFScUZa288l+wGuJEbXPG8xzfOcpM/Z3udc1DZ+kiZOiQdpKj4H1jKz1dL2rUcMw0vx\nFe/CXKBf8qWbuhR7zNHM9ALQx2JOX3WPv25qTls1j19b62a09XnadlXP5bKM5/ITd38IwN0fdPdd\niC9uJ4ZD4u4fufthxJDA0cC4LM8xsz/lxDC+lXUvcfb6EOLsT3WvWeZzhnhN1s/o1+dVHJfNAUQ9\ng7rI+ronPOP4f2Ycv5q753qWvqrXAIBkDuljwFXAz9x9DeBZYrhjfcj5eVT3uUrOhFzs7t2AnYgz\nNKmzfz8QgU1KVcFSZh/uz+jDKu5+ef08XRGRlaK4RHFJql+1jUv+p53kNehI9Z+fucC0jNd7VXc/\nqQ6Pndn2V0QiavO0tjt4FKbE3Re6+xnuviHQHzg9VcuhuveZuv3mr0z8JE2AEg/SJLj7XGIe3SiL\nAnZbEVnzB7Lc7RbgMjNbH8DMOpnZftUcez/xpfeYRfGeVmbW0czONbO9gNeIs9NnWxQi6gHsS8zB\nrKuTzeyXFsWSziOGPVbnduBEM9vBwipmtreZrWZmm5pZr+Q/pouJH4+K5DkfaWadkjMiC5K2Kqp+\nCI40s25m9hNiDuM4r/0yjl8Q8/LSPUnMYfwDMYc/2307mlmHtH0PAecn793axNy/bO85EBWazayL\nmV1PzCu8OPensIJqX/dqjn8dWGhRMKl90o8tzOw3OT5eVa9fSltiSOM8YJmZ9QN6V3Psysr6PKr7\nXJlZTzPb0mIo7HdEkJj6vL0NDEz+fkqAg2vowwPAvmbWJ3n8YoviVL+s92crIlJLiksUl9QmLqnC\nQ8DRZrZ18jqNBF5z90+qOf4ZYBMz+13yfrcxs9+Y2WZ1eOwVJO/F7cDVZvZTADNbx8z6JNf3MbOu\nZmZEbYjlxG9+te8z8Zu/l5mtZWY/J0ZGZLOy8ZM0AUo8SFNyGFG45nNiruJFaUPWqnItUSjneTNb\nSBQD2qGqAz3m1+0JfEAU8/mO+BJcm/ghWEr8oPcjMsM3AYM8VgeoqweJOfpzgI+JQk9VcvcyYunD\nG4BviCFvg5Ob2wGXJ/36L3EWYVhyW19gppl9T7weA5MholW5nygG9F+i0FBdqgkPB+61GCY3IOn7\nIuJMfRfg8SzP8QPih3hOcv/OxGtSBvwdeBd4kyyvE8mKFMT7N5VYfuw37v5uHZ5LTa97VccvJ87y\nb00U8fqKmPPbobr7ZBhFBDQLzOzMjLYXEu/J2KQvhxOf73qXw/Oo7nP1c2Ac8frPIuZCp+bOXkAU\nxvyGSAQ9WEMf5hJFv84lki1zgbPQ75aINB6KSxSX1BSXVNf2C8Tv4mNE/YWNyFLPIIkBeifHfE68\nJqlin/VhKPEevmqxCssLRPFIgI2T7e+J4qA3ufsUsr/P9xO1JT4hPlPZklj1ET9JE2Du1Y7sFZE8\nMbNPiGrD2QKUZsPMLgQ2cfcjC90XERERWZHiEhHJt5wKmomI1FUyZHMItV+NQkRERKReKS4RKQwN\nWRWRvDGz44gh8n9z9xcL3R8RERFpuRSXiBSOplqIiIiIiIiISN5oxIOIiIiIiIiI5I0SDyIiIiIi\nIiKSNyouWQtrr722b7DBBoXuhoiISF7MmDHjK3fvVOh+NFeKI0REpDnLFkco8VALG2ywAWVlZYXu\nhoiISF6Y2b8K3YfmTHGEiIg0Z9niCE21EBEREREREZG8UeJBRERERERERPImr4kHM+trZh+a2Wwz\nO6eK29uZ2SPJ7a+Z2QZptw1L9n9oZn1qatPMuiRtzE7abJt22wAze9/MZprZg2n7l5vZ28llfD5e\nAxEREakbxREiIiLNQ94SD2ZWBNwI9AO6AYeZWbeMw4YA37h7V+BqYHRy327AQGBzoC9wk5kV1dDm\naODqpK1vkrYxs42BYcDO7r458Me0x1/k7lsnl/71+wqIiIhIXSmOEBERaT7yOeJhe2C2u89x96XA\nw8B+GcfsB9ybXB8H7GFmlux/2N2XuPs/gdlJe1W2mdynV9IGSZv7J9ePA250928A3P3LPDxXERER\nqV+KI0RERJqJfCYe1gHmpm1/muyr8hh3XwZ8C3TMct/q9ncEFiRtZD7WJsAmZvaymb1qZn3T7l9s\nZmXJ/v2pgpkdnxxTNm/evFyet4iIiKw8xREiIiLNREtYTrM1sDHQA/gl8KKZbenuC4D13f0zM9sQ\nmGxm77r7x+l3dvfbgNsASkpKvGG7LiIiIgWmOEJERGQl5XPEw2fAumnbv0z2VXmMmbUGOgDzs9y3\nuv3zgTWSNjIf61NgvLuXJ8Mt/0EEELj7Z8m/c4CpwDZ1e6oiIiJSzxRHiIiINBP5TDy8AWycVIlu\nSxR5yqz4PB44Krl+MDDZ3T3ZPzCpVt2F+IF/vbo2k/tMSdogafOp5PqTxFkKzGxtYsjkHDNb08za\npe3fGXi/Pl8AERGRfCkthVGj4t9mSnFETVrAh0BERJqHvE21cPdlZnYKMAEoAu5y95lmdglQ5u7j\ngTuB+81sNvA1EQCQHDeW+AFfBpzs7ssBqmozecihwMNmdinwVtI2ybG9zex9YDlwlrvPN7OdgFvN\nrIJIwFzu7ko8iIhIo/fKK9CzJ5SXQ3ExTJoE3bsXulf1S3FEDSZPht69wR3atWueHwIREWk2LJL8\nkouSkhIvKysrdDdERKQFmzIFBg2CTz+N7aIiGDEChg1b+bbNbIa7l6x8S1KVeo0jjj8ebr89rpvB\n0KEx+kFERKRAssUR+ZxqISIiIvVkxow4wd2rFyxdCm3aRNKhbVvo0aPQvZMGd/TRMdzFLEY9jBkD\nZ58NX39d6J6JiIj8j5awqoWIiEiTVFoKjz4K77wTI+s7doz/X550Erz1FkydGkkHjbBvgbp3jw/F\n1KmwySYwfjxcdRXcdhucdVbc/tpr+oCIiEijoKkWtaCpFiIi0lCeegoOOgiWL4/tY46Bq6+G1VfP\n32NqqkV+5T2OeO89OP/8+PBAjIZorkVARESk0dFUCxERkSZi/vw4YZ2edCgqgq5d85t0kGZgiy3g\nySfhxBNj2x0WLYJrroGKisL2TUREWjQlHkRERBqBH36Ayy6DDTeEP/856jkUF6uOg9TBoEHQvj20\nahWjHsaOhW22gWefjWSEiIhIA1ONBxERkQJ68cVINEyfDt98A/37RwJiiy2ixoPqOEitde8e0yum\nToXddoN//zumYOy9N+y6Kxx5ZAyt0QdLREQaiGo81IJqPIiISH2pqIhlMC++OE5Ct2oFN98cqyQW\nimo85FdB44ilS+HOOyMBkVr5org4ClQq+SAiIvVANR5EREQaCXeYMAFKSmD48MqR72ZxElokL9q2\njeVQTjstPmwAixfDKafEiAgREZE8UuJBRESkgbz+OuyxB/TtG9MqLrwwpuKrjoM0mPTiIa1bw7vv\nxnKcZ5yhzJeIiOSNajyIiIjkUWkpPPoovPVWTLnv1Amuuy6mVLRrF0kI1XGQBpNe/6FHD1hnnRh6\nc801cMcdMHAgdO4cCQp9IEVEpJ6oxkMtqMaDiIjUxvjxcOCBlctiDhkCV18Nq61W2H5VRzUe8qtR\nxxEzZ8ZUjOnTY7tNm0hQ7LprYfslIiJNhmo8iIiINKAFC+Dcc+GggyqTDkVFsNFGjTfpIC3c5ptD\nv35R5RSgvBwOOCCG6+gklYiIrCQlHkREROrJ4sUwZkwkGEaNgp49K6fTq4aDNHo9esT8n9QHdvXV\nYcAA2GGHWP1CRESkjlTjQUREZCW99BJce22MUv/iC+jTJxIP22wTNR5Uw0GahMz6D9tvDw88ABdc\nEFVR+/SBww6Dzz/XB1pERGpFNR5qoVHPzRQRkQbnHiMczjorrptFAuLUUwvds7pRjYf8arJxxOLF\ncOONcPHFsHBh7CsujlEQSj6IiEhCNR5ERETq2euvx1SKM8+snALfqhV8/31h+yVS74qLY7nNP/4x\nsmsQyYgzz4Svvips30REpElQ4kFERKQWPvoIDjkkpr3PmhX/H2vfXnUcpAXo16+yaElRUcwj2mgj\nGDkSfvyx0L0TEZFGTDUeREREcvDXv8KIEVBWFv/3uuiiSDqstlqsXqE6DtLsZdaA6NABhg2D886L\nqRiDBsGqq0KvXvpDEBGRFSjxICIiksXChTHC/K67Yrt1a3jkEdh778pjunfX/7Okhcj8sD/1VFRX\nPekkuPzy2Ne2LUyZAjvtVJg+iohIo6OpFiIiIlUoL4ebb4auXSuTDhD1HP7+98L1S6TR2WWXWO2i\nVRJWLl0Khx8Or75a2H6JiEijocSDiIhIGnd4/HHYYgv4/e/hV7+C229XHQeRrHr2hHbt4o+kTZsY\nKtS9Oxx8MPzjH4XunYiIFJimWoiIiBB18u65B155Bd57D7p1g6efjikVZrD55qrjIFKtzPoPW24Z\na81eeSU8+STst19k8fbZR39AIiItkHlqDTCpUZNdf1tERLJ66CE48kioqIjtYcPgkkuinkNLkm39\nbVl5LTKO+OILOPlkeOyx2G7dGp57DvbYo7D9EhGRepctjtBUCxERabH++9+oiXfEEZVJh6KiWKmi\npSUdRPLiZz+D7barrP+wbBnsvz/cdltcFxGRFkGJBxERaXG+/x4uvjgKR95xRyyHWVysGg4iedGj\nR2X9h3btoEsXOOGEmI4xfnwUVhERkWZN53NERKTFmD49pp1PmwbffBN170aNigREaalqOIjkRWb9\nhx13jGU4hw6N2g+77gqDBsG8efoDFBFpplTjoRZa5NxMEZFmwD1q3J1zTlxv1QpuuQWOO67QPWtc\nVOMhvxRHZCgvhzvvhHPPjUwgxNCjyZOVfBARaYJU40FERFqsN96Ik6hDh1aO6DaDr74qaLdEpE0b\nOPFEOO20+KMEWLwYzjwT5s8vbN9ERKReKfEgIiLN0j//CYcdBttvD7NmwRlnQPv2quMg0uj06VNZ\nZKWoKOY9bbRRDFNavLjQvRMRkXqgGg8iItKsTJgAl14Kr74aJ1TPOw/OPhtWXz2KSKqOg0gjk1kD\nYrXVYojS2WfDjTfC4MGRLezZU3+4IiJNlBIPIiLSLCxZAmedBddfH9tFRfDww9C/f+Ux3bvr/y0i\njVLmH+df/xq1Hk46KZaggVgRY8oU/RGLiDRBmmohIiJNmjuMHQubbVaZdEiZObMwfRKRetCrFxx1\nVFSDhcguHnNMzJ0SEZEmRYkHERFpsl56KU5+HnpojM6++mrVcRBpVnr2jJEORUUxd+rf/4Ytt4Tf\n/x6+/LLQvRMRkRxpqoWIiDQppaUwbhzMmAHTpkHnznDXXTBoUPzfZIcdVMdBpNnIrP/QtSuMGAE3\n3wwPPACHHw6/+AX07q0/eBGRRsw8tbaY1Ejrb4uIFNZzz8G++8KyZbF93HExymGVVQrbr+Yi2/rb\nsvIUR9Sjf/wjvgBefDG227aNmhA771zYfomItGDZ4ghNtRARkUZvyRL485/hgAMqkw5FRdCli5IO\nIi3SJptA376V9R+WLoVDDonikyIi0ugo8SAiIo1WeuHIM8+ErbeG4mLVcBAR4gsgVf+hbVuoqIiC\nlP37wwcfFLp3IiKSRjUeRESk0Skthfvug+nTY2WKrbaC55+H3/42blMNBxH5n/oP22wD114LI0fC\nFlvACSfA3nvDO+/oC0NEpMBU46EWNDdTRCT/xo2DgQNh+fLYPvdcuOSSOKkp+aUaD/mlOKKBzJsH\nF18cBSgrKsAshkpNmqTkg4hIHhWsxoOZ9TWzD81stpmdU8Xt7czskeT218xsg7TbhiX7PzSzPjW1\naWZdkjZmJ222TbttgJm9b2YzzezBtP1HmdlHyeWofLwGIiKSmwUL4KyzVkw6FBXBqqsq6dBSKY6Q\nOunUCW64Af7wh9h2h0WLYjSETriJiBRE3hIPZlYE3Aj0A7oBh5lZt4zDhgDfuHtX4GpgdHLfbsBA\nYHOgL3CTmRXV0OZo4OqkrW+StjGzjYFhwM7uvjnwx2T/WsBFwA7A9sBFZrZmvb8QIiKSVXl5/B+h\na9coINm3r+o4iOIIqQeHHALt20cBSjN45BHYaaeYryUiIg0qnyMetgdmu/scd18KPAzsl3HMfsC9\nyfVxwB5mZsn+h919ibv/E5idtFdlm8l9eiVtkLS5f3L9OOBGd/8GwN2/TPb3ASa6+9fJbROJ4ERE\nRBrAK6/AoEGRcDj11Kjj8Oab8MwzsSreiBEaGd3CKY6QlZOqAXHppbHs5l13wb/+FcmHQw+Fxx6D\nUaOUiBARaQD5LC65DjA3bftT4qxAlce4+zIz+xbomOx/NeO+6yTXq2qzI7DA3ZdVcfwmAGb2MlAE\nDHf356rp3zqIiEje3XsvHHNM5fTrK66IVSvM4vbu3ZVwEMURUg/Sv0x22QUGDIArr4TLL48lc1T/\nQUSkQbSE5TRbAxsDPYDDgNvNbI1c72xmx5tZmZmVzZs3L09dFBFpGf7zHxgyBAYPjqQDxCjoZcsq\nkw4ijYziiOZklVVg+HD4059iO1X/4fLL44tIRETyIp+Jh8+AddO2f5nsq/IYM2sNdADmZ7lvdfvn\nA2skbWQ+1qfAeHcvT4Zb/oMIIHLpH+5+m7uXuHtJp06dcnjaIiKSadGiGO288cZw//1RQLJ9e9Vx\nkKwUR0j+9O9f+SXUqhWMHx/zvf72t0L3TESkWcpn4uENYOOkSnRbosjT+IxjxgOpKtAHA5M91vcc\nDwxMqlV3IX7gX6+uzeQ+U5I2SNp8Krn+JHGWAjNbmxgyOQeYAPQ2szWTYlC9k30iIlJPXn45plJv\nsAFccAH06QPvvw8PPRQjm1XHQbJQHCH5k6r/MGIETJ8OTzwRlW732isq3D7wgOo/iIjUo7zVeEjm\nWp5C/AgXAXe5+0wzuwQoc/fxwJ3A/WY2G/iaCABIjhsLvA8sA0529+UAVbWZPORQ4GEzuxR4K2kb\nKgOD94HlwFnuPj9pawQRhABc4u5f5+v1EBFpaW67DU46qbKOw403wu9/X3m76jhINoojJO8yv4T2\n2gtuuimypBMmxBdXu3ZR7VZfViIiK8Vc6xnnrKSkxMvKygrdDRGRRu1f/4KhQ2PlupSiojixOGxY\n4folNTOzGe5eUuh+NFeKI5qICy6Ayy6L+g8QIyCefDKSECIiUq1scURLKC4pIiINYOFCOPdc2HTT\nmC599NGq4yAiTdBee8VKF6n6D889B926weOPVyYjRESkVvK5nKaIiDRzpaUxCvmHH+Cuu+CLL+DI\nI2HkSFh3XTjuOJg6NZIOGqksIk1Cqv5D6str4UI4/XQ46CDYffdYluc//9EXm4hILSjxICIidVJa\nCj17wpIlsb3FFjHSYfvtK49RHQcRaZIyv7zefhvuuAPOOSeGc6n+g4hIrWiqhYiI1NrHH8Oxx1Ym\nHVq1gsMOWzHpICLSbLRuDSeeCKedFkkHd1i8GIYPj39FRCQrJR5ERCRn330XhSO7dYM5c6BNm5gG\n3a5djH4QEWnW+vVbsf7D88/Dr34Fjz6q+g8iIllknWphZr8klqbaFegMLALeA/4K/M3dK/LeQxER\nKahUHYeFC+Huu+HLL2OK82WXxQoWquEg1VEcIc1OZv2HRYvgT3+CAQNgl13gmmtg6VJ9MYqIZKg2\n8WBmdwPrAM8Ao4EvgWJgE6AvcJ6ZnePuLzZER0VEpOFl1nHYaiv461+hJFkoqXNnxdVSNcUR0mxl\n1n94882ornveefHlWFQU+9u2jSSFviRFRLKOePizu79Xxf73gMfNrC2wXn66JSIihTZnTqxKkV7H\n4dBDK5MOIjVQHCEtQ1FRfFkOGAD9+8OLSS5t8WJ44QUlHkREyFLjoZpgIf32pe4+u/67JCIihbRw\nIQwbBpttFkUkVcdB6kJxhCEIu1wAACAASURBVLQ4HTrA5ZdHDQiImg833wyPPab6DyLS4tWpuKSZ\n/a2+OyIiIoVVUQH33AObbBKx86GHwuzZMG0ajBihEcNSfxRHSLPVvXsUxRk5Eq67Djp2hIMPjqzt\n228XunciIgWTrcbDttXdBGydn+6IiEgh3HprJBc++wx23BGefBJ22CFuW2cdJRyk9hRHSIuVXgPi\npJPg9tvhggtg221h331hiy1gn330xSoiLUq2Gg9vANOIACHTGvnpjoiINKS5c2HIEJg4MbbbtoWr\nrqpMOoisBMURIq1bR/Jh4EA48UQYOxbGj4crr4z6D7vtVugeiog0iGyJh1nACe7+UeYNZjY3f10S\nEZF8W7Qo4t7LL4fycjCLKcjLl0ddtJ13LnQPpRlQHCGSsuaasPXWMG5czGsrL4eDDoJ774W99ip0\n70RE8i5bjYfhWW4/tf67IiIi+eYOjz4Kv/oVXHRRjPZ95JGohVZUFCMeevQodC+lmRiO4giRSj16\nRJXe1Jdt+/aw995x+fDDQvdORCSvqh3x4O7jstz2ZH66IyIi+VBaCg8+CC+9FPXNfv1ruO8+2H33\nuH3SJJg6NeJiTTuW+qA4QiRD9+4rftlutx1cfz1ccknUfTj4YNh4Y+jXT1/EItLsmNdieR8ze8bd\n98ljfxq1kpISLysrK3Q3RERq5dlnY2n55ctj++yzo+B6UVFh+yWNj5nNcPeSPLavOEJxhGT64gs4\n7jh4+unYbtMmVsbYZZfC9ktEpJayxRG1XU5znXroj4iINIDycrj22phGnEo6FBXBGmso6SAFozhC\nJNPPfhYjHFolYXl5eaxn/PLLhe2XiEg9qm3i4a289EJEROrVxIlRx+yPf4SttlINB2k0FEeIVCWz\n/sPSpTHi4Ygj4NNPC907EZGVlm1Vi//h7sfkqyMiIrJySkvh8cfhtddg+nTYcEN48smYZvHqq6rh\nIIWnOEKkGpn1H7baKpYduvLK+CI/99xYbqi0VF/kItIk1Vjjwcx+6+4TG6g/jZrmZopIYzVpEvTt\nC8uWxfaJJ8I118QJNJFc5aPGg+KISoojpNb++U8488zIKpvFpV27+NJX8kFEGpk613gws0HA6Xnp\nlYiIrDR3eOghOPDAyqRDURGst56SDlJ4iiNEVlKXLvDYY3DMMfGFX1EBixbB2LGF7pmISK1Um3gw\ns/OAo4GDGq47IiKSq7ffht12g8MPh5//fMXpwarjIIWmOEKkHh17LLRvX1mA8vrr4fTT4dtvC9sv\nEZEcZavxcCGwmbv/2FCdERGR7EpLY3nMmTPhqadgrbXg9tvh6KPh9ddVx0EaFcURIvUlvQbEr38d\ndR+uuQb+8peoBbHJJvDii/oBEJFGK1vi4WhgnJn1dvevGqpDIiJStenTYY89YqU1gEMOgVtvhTXX\njO3u3RVvSqOiOEKkPqV/ye+1F5xwApx6akzDUP0HEWnkqp1q4e4PAmcCTzZcd0REpCovvggDBlQm\nHYqKYJttKpMOIo2N4giRPNtuO3j55chCp+o/LF4MzzxT6J6JiPyPrMUl3X0ycGID9UVERDJ89lnU\ncNh994gr27ZVHQdpOhRHiOSZGfzpT1H/wSx+KK65Ji6pTLWISCOQNfEA4O7vNURHRESk0pIlMW13\n001jFbULLoA5c2J674gRGkkrTYfiCJE8S9V/uOyyWOZot90iGbH11jB5cqF7JyICZKnxYGZHuvsD\nZlblMljuPiZ/3RIRaZlKS6Nuw6RJ8OmnsN9+MGYMbLhh3K46DtJUKI4QaUDpPw6HHgpPPw1//GMU\nBurZE37zG9h/f/2AiEjBZCsuuUry72oN0RERkZZu3LiIFysqYsTsmDFx0kqkiVIcIVIIZtC/P/Tu\nDaedFksfTZkSPyoTJ2qenogURLWJB3e/Nfn34obrjohIy/PjjzB6dIySraiIfa1aRY0wkaZKcYRI\ngRUXQ5cuURho+XJYtgwOPBDuuw/22afQvRORFqbGGg9mdoWZrW5mbcxskpnNM7MjG6JzIiLNmTs8\n8QR06waXXAK9ekWcqOKR0pwojhApoB49KqsSt2sHq68O++4bl48/LnTvRKQFqTHxAPR29++AfYBP\ngK7AWfnslIhIc1ZaCqefDjvuGCefVl8dpk2D55+POmAqHinNjOIIkUJJFZ4cMSKmW/zjH3DVVVGp\nuFu3qFw8ZQqMGhU/TiIieZKtxkPmMXsDj7r7t2aWxy6JiDRfkyZB374x4hWi9teVV0Lr5JtWxSOl\nGVIcIVJImT8sZ5wR6zSffTZcemnM8zOLERHKeotInuQy4uEZM/sA2A6YZGadAM08FhGpBXd45JEY\n4ZBKOhQVwU9/Wpl0EGmmFEeINDa/+AXcfz8cf3z8QFVUwKJF8Oijhe6ZiDRTNSYe3P0cYCegxN3L\ngR+A/fLdMRGR5uL992HPPWHgQPj5z+Okkuo4SEuhOEKkERs8GNq3j4rGANddB+ecA99/X9BuiUjz\nU+N5NjMrBgYDu5iZAy8BN+e5XyIiTVppKUyYENNpH30UVl0VbrwRTjgBXn89ptf26KERrdL8KY4Q\nacRSNSCmToWttoofrNGj4S9/ieU311knihDpB0tEVpK5e/YDzMYCC4EHkl2HA2u4+yF57lujU1JS\n4mVlZYXuhog0cq+8Aj17wtKlsb3vvnDnndCpU2H7JVITM5vh7iX13KbiiITiCGkSXnkFTj4Z3n67\nciSE6j+ISA6yxRG5zCzewt27pW1PMbP366drIiLNy8yZcOSRlUmHoqKI05R0kBZMcYRIU7LTTlBW\nBgccAE8/HfsWL45hfEo8iEgd5VJc8k0z2zG1YWY7AErXi4ikWbgQzjoLtt4avvoK2rRRHQeRhOII\nkaamqAiGDYPi4ljxwh1uugnGjYvrIiK1lEviYTvgFTP7xMw+AUqB35jZu2b292x3NLO+Zvahmc02\ns3OquL2dmT2S3P6amW2QdtuwZP+HZtanpjbNrEvSxuykzbbJ/sFmNs/M3k4ux6bdZ3na/vE5vBYi\nIitIrVbxq1/F0uhHHQVz5sSU2BEjNDJVBMURIk1T9+4weXIst3nrrdC5MxxyCPTpAx9+WOjeiUgT\nk8tUi751adjMioAbgd8CnwJvmNl4d08fXjkE+Mbdu5rZQGA0cKiZdQMGApsDnYEXzGyT5D7VtTka\nuNrdHzazW5K2U8WrHnH3U6ro5iJ337ouz09EWrbS0qjBNX16jEjdZht47DHYMTmvu/baSjiIJBRH\niDRV3btX/pgdcwzcfDOcfz5suWUs1bTRRtC7t37wRKRG1SYezGxVd//e3f+V7ZgsbW8PzHb3Ocmx\nDxPLZ6UHDPsBw5Pr44AbzMyS/Q+7+xLgn2Y2O2mPqto0s1lAL6JgFcC9Sbuqmi0i9W7y5Djhs2xZ\nbJ9xRhQBLyoqbL9EGhPFESLNTOvWcOqpMGAAHH003H9/7B85EqZMidoQIiLVyDbV4ikz+7OZ7WZm\nq6R2mtmGZjbEzCaQ/SzGOsDctO1Pk31VHuPuy4BvgY5Z7lvd/o7AgqSNqh7rIDP7u5mNM7N10/YX\nm1mZmb1qZvtneS4iIrjDE0/AgQdWJh2KiqBjRyUdRKqgOEKkOfrZz2DXXStXvFi6NBIRH39c2H6J\nSKNWbeLB3fcAJgEnADPN7Fszm08sh/Vz4Ch3H9cw3VwpTwMbuPtWwETiLEbK+slyH4cD15jZRpl3\nNrPjk6CibN68eQ3TYxFpdD7+GPbeO5IOHTvGymIqHilSPcURQXGENEs9elT+ELZpA3Pnwuabw8UX\nxwoYIiIZstZ4cPdngWfr2PZnQPpZgV8m+6o65lMzaw10AObXcN+q9s8H1jCz1snZiv8/3t3npx1/\nB3BFasPdU8fMMbOpwDbACulad78NuA1i/e0cnreINBOlpfDCC/Cvf8EDD0RsNWYMnHJK1HWYOjVi\nL01tFama4gjFEdJMde8e1ZNTP4TrrRfzDocPjykYJ58cCQj9SIpIIpfiknX1BrCxmXUhfrwHUjl3\nMmU8cBRR4fpgYLK7e1IZ+kEzG0MUhdoYeB2wqtpM7jMlaePhpM2nAMzsF+7+n+Tx+gOzkv1rAj+6\n+xIzWxvYmbRgQkRattJS6NkTliyJ7T33hHvvjaLesGK9LRHJC8URIo1Z5g/hww/DscfCkCFw+umx\nr7g4CiPpB1OkxctlOc06Sc4YnAJMIH6kx7r7TDO7xMz6J4fdCXRMij6dDpyT3HcmMJYoIPUccLK7\nL6+uzaStocDpSVsdk7YBTjOzmWb2DnAaMDjZvxlQluyfAlyeUSlbRFqouXPhuOMqkw6tWkGvXpVJ\nBxHJP8URIk3QnntG8sEsthcvjrWly8sL2y8RKThz16i/XJWUlHhZWVmhuyEieVJeDtdeGyNFy8uj\nmGRFRdRxmDRJJ2yk+TOzGUnNAskDxRHSIpSWwh57RNHJ1A/p5pvHUpy77lro3olIHmWLI2oc8ZBU\npN68/rslItJ4vPQSbLstnHVWTLH44AOYNi1O1CjpIFJ3iiNEWphU/YcRI+LH9amn4PvvYbfdYPBg\n+PLLQvdQRAoglxoPs4DbkqJNdwMPufu3+e2WiEjDePZZuOACePPNqI315JPQv3+MEu3SRQkHkXqg\nOEKkpcms/7DHHnDZZXDVVZGIOO446NAh5jHqh1akRahxxIO73+HuOwODgA2Av5vZg2bWM9+dExHJ\nl4oKGDo0lsh8801o3Rruvhv2269yaqqIrDzFESLCKqvAyJHwzjuw4YZw5ZVw/vkxxLC0tNC9E5EG\nkFNxSTMrAn6VXL4C3iEKMD2cx76JiOTFO+/ALrvAFWn1593htdcK1yeR5kxxhIgAsNlmcNBBUbUZ\noorzmWfCd98Vtl8ikne51Hi4GvgA2AsY6e7buftod9+XWK9aRKRJWLgwlhnfbjv46KOYYtG+PRQV\nRQHJHj0K3UOR5kdxhIisoGdPaNcufnxbt4ZXXomExKOPxlkAEWmWcqnx8HfgfHf/oYrbtq/n/oiI\n1LtXXoEbb4Tnn4evvoLjj4dRo2CttaBfP5g6NZIOmmYqkheKI0SkUqr4ZOrHt6gITjwRBgyAPn1g\nyBCYPVs/zCLNTC6JhyPd/e70HWY2yd33UHEoEWnsHnssYpmKiqjdcNttUdMqJbP+lYjUO8URIrKi\nzB/f11+Hm26CYcNgwoT4wS4u1rJSIs1ItVMtzKzYzNYC1jazNc1sreSyAbBOQ3VQRKQuli6Fyy+H\ngQMj6QAxpfSrrwrbL5GWQnGEiOSsdWs47bS4QEy5WLQI7rqrsP0SkXqTrcbDCcAMohDUm8n1GcBT\nwA3575qISN1Mnw7bbhsnTnbeOU6aqI6DSINTHCEitdO/fxRfatUqRj3ccQccdRR8+WWheyYiK6na\nqRbufi1wrZmd6u7XN2CfRETqZP58OPvsOEGy/vrw9NOwzz6xUpfqOIg0LMURIlJr6fUfUtdHj44f\n9NGjo/5Dq5wW5RORRsa8muqxZtbL3Seb2YFV3e7uj+e1Z41QSUmJl5WVFbobIpLhlVdgzBiYOBF+\n/DFWrrjgglg2XERyZ2Yz3L2kntpSHJFBcYRIHcyaBSedBNOmwZZbxqoYAwfqTIJII5QtjshWXHJ3\nYDKwbxW3OdDiAgYRaXweegiOPDLqOLRqBffeG9siUnCKI0Rk5W22GUyZAhdeCJdeCu++G4UoJ0yA\nXr0K3TsRyVG2qRYXJf8e3XDdERHJzeLFMHIkXHZZZfFIM5g7t7D9EpGgOEJE6o0Z/OQnUbBp+XJY\ntgwOOggefDDWxRaRRq/axIOZnZ7tju4+pv67IyJSs0mTYtTlRx/Fkt/TpkF5uYpHijQmiiNEpF71\n6BE/9EuXxioYHTrAXnvBIYfANddA586F7qGIZJFtqsVqDdYLEZEalJbCM8/AjBkxurJr16jpsOee\nKh4p0kgpjhCR+pNeeLJHj1i+6sorY/rFhAkwahRstVUsbaWAQKTRqba4pPwvFYUSKYyXX45aUuXl\nsT14cEzvbN++oN0SaXbqs7ik/C/FESJ5MHt2DIN84YXKFS/atYskhZIPIg2qTsUlzexsd7/CzK4n\nikCtwN1Pq8c+iohUadasKBaZSjoUFcEmmyjpINLYKY4QkQbRtSs8/3ysdDF2bOxbvDj2KfEg0mhk\nm2oxK/lXqXkRaXCLF8eoyVGjoLgY2rSJIpKq4yDSZCiOEJGGYQZ//CM8/XQEEO5wyy2www7Qt2+h\neyciZF/V4unk33sBzGz12PSFDdQ3EWmhpkyBE0+Ef/wDjjgCxoyBjz9WHQeRpkRxhIg0qPQaEB06\nwPXXx4oXhx4axSd//vNC91CkRcs24gEAMysB7iaKRJmZLQCOcfcZ+e6ciLQszz0H550Hb74JG24Y\ntaJ6947bfvpTJRxEmiLFESLSYLp3rwwWhgyB0aNj3e0JE+D44yMh0bOnAgqRAqgx8QDcBfze3acD\nmNkuRACxVT47JiIthztcdBGMGBHbrVvDHXdEbCAiTZ7iCBFpeO3awYUXxoiHww6DK66I/cXFMHmy\nkg8iDaxVDscsTwULAO7+ErAsf10SkZbk449jVEMq6QCRiHj11cL1SUTqleIIESmcTTeFgw+OOhAQ\nNSDOOy/+FZEGU23iwcy2NbNtgWlmdquZ9TCz3c3sJmBqg/VQRJql8vIYAbnFFvDaa3DGGbFSRVGR\nCkiKNAeKI0Sk0ejZM0Y6FBXFZcoU+PWvYdq0QvdMpMXINtXizxnbF6Vd/59lsUREalJaGjWfOnWK\nmk9//zsceCBcdx2ssw4cdJAKSIo0I4ojRKRxSC882aMHLFwYVax79IhaEIccEgWmFICI5E22VS00\nu1pE6k1pKeyxR+UqV2uvDU88AfvvX3lMek0oEWnaFEeISKOSGWS89x5cfDFcdRXceSe0ahV1ISZN\nUjAikgfVJh7M7Eh3f8DMTq/qdncfk79uiUhzc/PNsGhRXDeD3/9+xaSDiDQviiNEpFH7yU9izueP\nP8INN0BFRQQqTzyhxINIHmQrLrlK8u9q1VxERGr03//CgAFw//2RcGjVKqZZ9u1b6J6JSJ4pjhCR\nxu/ww6PIVKvkv0XXXw/XXAPLlxe2XyLNTLapFrcm/17ccN0RkebCHe66C848M04mjBgBu+0GL7+s\nKZQiLYHiCBFpEtLrP2y6aUy7+NOf4MEH4fbbowiliKy0bMUlATCzK4BLgUXAc8S6239y9wfy3DcR\naYJKS2HcuCgY/dZbkWy47bb4LYfYFpGWQ3GEiDR66fUfDjgAHnkE/vAHKCmJMyh9+kSAozMnInWW\nbapFSm93/w7YB/gE6Aqclc9OiUjTNH067L47jBkTSYehQyMBkUo6iEiLpDhCRJoOMxg4EGbNgt/9\nDi6/HHr1gvPPjyrZpaWF7qFIk5RL4iE1KmJv4FF3/zaP/RGRJuqNN6KWQ3l5bBcVQYcOlVMmRaTF\nUhwhIk3PWmvFnNEhQ2L+aEVFLM31t78VumciTVIu/yV4xsw+ALYDJplZJ2BxfrslIk3FDz/AGWfA\njjvCsmXQtm0kHdq2jRGJItLiKY4QkaZryJAoPmkWCYgbb4THHit0r0SanBoTD+5+DrATUOLu5cAP\nwH757piINH7PPw9bbBFTK44/HmbPjtpMI0ZoGWwRCYojRKRJSxWfvOwyuPtuWH99OPhgOPBA+Pzz\nQvdOpMnIpbjkIcBz7r7czM4HtiWKRP03350TkcbpuefgvPPgzTdhk03gxRdh113jtvT6TCIiiiNE\npMlLD26OPBKuvhouvBA22wxOPDHmlvbsqQBIJItcplpc4O4LzWwXYE/gTuDm/HZLRBoj9xjN0K9f\nJB1at4Zbb61MOoiIVEFxhIg0H61bw1lnwbvvwsYbwxVXxNmYXr1UeFIki1wSD8uTf/cGbnP3vwJt\n89clEWmMPv0U9tsvEvwp7vqNFZEaKY4Qkeana9eYbmEW24sXw8UXV1bZFpEV5JJ4+MzMbgUOBZ41\ns3Y53k9EmoGKCrjlFujWDV54AU49NWosqYCkiORIcYSINE89e0JxcQRFRUUwYQLssEMMCxWRFdRY\n4wEYAPQFrnL3BWb2C7T+tkizVloaRSI33BBuuilqOPTqBbfdBhttBIcdFrf36KHpjCJSI8URItI8\npQpPpoKi//wHTj4Ztt8ezjwT+vaNoEoBk0jNiQd3/9HMvgR2AT4CliX/ikgzVFoKe+wRIwbdYdVV\n4c474eijK0cTqoCkiORKcYSINGuZQVHPnpF0GD066j+YQbt2Wu5LWrwahzqa2UXAUGBYsqsN8EA+\nOyUihfPgg7BoUSQdIKZWHHNMZdJBRKQ2FEeISIuy5ppxxuaYYyKYqqiIsznPPVfonokUVC5zLA8A\n+hPrbuPunwOr5dK4mfU1sw/NbLaZnVPF7e3M7JHk9tfMbIO024Yl+z80sz41tWlmXZI2Zidttk32\nDzazeWb2dnI5Nu0+R5nZR8nlqFyek0hztXhxFGW+6abYbtUqajnsu29h+yUiTZ7iCBFpeY49NgIp\ns0hA3HwzPPtsoXslUjC5JB6WursDDmBmq+TSsJkVATcC/YBuwGFm1i3jsCHAN+7eFbgaGJ3ctxsw\nENicmBd6k5kV1dDmaODqpK1vkrZTHnH3rZPLHcljrAVcBOwAbA9cZGZr5vLcRJqbl1+GrbeGkSNh\n0KBIyl96qUYFiki9UBwhIi1Pqv7DZZfB7bfDT38Ke+8NRx4JX31V6N6JNLhcikuOTapRr2FmxwHH\nALfncL/tgdnuPgfAzB4G9gPeTztmP2B4cn0ccIOZWbL/YXdfAvzTzGYn7VFVm2Y2C+gFHJ4cc2/S\nbrZ1wvsAE93966StiURw8lAOz02kySstjeLLM2fCY4/BeuvFdu/ecXufPtnvLyKSI8URItIypdd/\n+N3vYNSoSEQ8/3zMZS0qipoQOssjLUAuxSWvMrPfAt8BmwIXuvvEHNpeB5ibtv0pcVagymPcfZmZ\nfQt0TPa/mnHfdZLrVbXZEVjg7suqOB7gIDPbDfgH8Cd3n1tN/9LvA4CZHQ8cD7DeeutleboiTUdp\nafzOLVkS2wcfDHffHYUkRUTqk+IIxREiQhSYHD4cDjoIBgyACy+M/cXFMHmykg/S7GWdapEMS5zi\n7hPd/Sx3PzPHYKExeRrYwN23AiYSZzFy5u63uXuJu5d06tQpLx0UaUgLFkSSPZV0aNUKtt1WSQcR\nqX+KIxRHiEiGLbeM6Rapqt2LF8NVV1VW9RZpprImHtx9OVBhZh3q0PZnwLpp279M9lV5jJm1BjoA\n87Pct7r984khnK0z9uPu85OhlgB3ANvVon8izcr48dCtG7z1FrRuHSP82rWL5aVFROqb4ggRkSr0\n6hUjHYqK4gzQ44/HHNdPPil0z0TyJpfikt8D75rZnWZ2XeqSw/3eADZOqkS3JYo8jc84ZjyQqgJ9\nMDA5KUA1HhiYVKvuAmwMvF5dm8l9piRtkLT5FICZ/SLt8foDs5LrE4DeZrZmUgyqd7JPpNmZNw8O\nPxz22w86dYLXX4cXX4QRI1RAUkTyTnGEiEi6VOHJESMiILvxxpgHu8UWcb2iotA9FKl3uRSXfDy5\n1Eoy1/IU4ke4CLjL3Wea2SVAmbuPB+4E7k+KPn1NBAAkx40lCkgtA05OzppQVZvJQw4FHjazS4G3\nkrYBTjOz/kk7XwODk8f42sxGEEEIwCWpAlEizcUrr8B118UqFT/+CJdcAkOHQtu2cbsSDiLSABRH\niIhkSi88ufPOseLFCSfAKafAI4/AySfDnDkxLFUBmzQD5jnMJ0rOCvyKWArrQ3dfmu+ONUYlJSVe\nVlZW6G6I5OTpp2H//SNpbgb33RdTCkVEqmNmM9y9JA/tKo5AcYSI1MAd7r03inF9/30EcMXFGp4q\nTUa2OKLGqRZmthfwMXAdcAMw28z61W8XRaS+pH6zBgyoHKnXqhXMnZv9fiIi+aA4QkQkR2YweHAk\nHiCCukWLYOzYgnZLpD7kUuNhDNDT3Xu4++5AT+Dq/HZLROri009jpN7gwbDxxpV1i9q2VQFJESkY\nxREiIrWx777Qvn2cOQK44QYYORLKywvbL5GVkEviYaG7z07bngMszFN/RKQO3OH222HzzWHaNLj2\nWnj77VgWWgUkRaTAFEeIiNRGqvjkpZfCM8/AgQfCeefBDjvAO+8UuncidZJLcckyM3sWGEvMzTwE\neMPMDgRw91oXjBKR+vP443D22fDxx9CzJ9xxB2y4YdyWXrdIRKRAFEeIiNRWehC3995w6KFw0klQ\nUgKDBkGXLrDHHgr0pMnIJfFQDHwB7J5szwPaA/sSAYQCBpECqKiIhMOf/xzbbdrE6IZU0kFEpJFQ\nHCEisrIOPBB23x2OOALuuiv2XXopTJmi5IM0CTUmHtz96IboiIjkbs4cGDIEpk6t3FdREUtB77xz\nwbolIvI/FEeIiNSTjh0j+TBxYgR+S5bAuefGuunt2hW6dyJZ5bKqRRczG2Nmj5vZ+NSlITonIiuq\nqIDrr4ctt4Q334Rhw6L2kApIikhjpThCRKQe9egRSYaiorhMnQrbbQdvvFHonolklctUiyeBO4Gn\ngYr8dkdEqjN7NhxzDEyfDv36wa23wrrrRuHjqVPjd0gj7USkEVIcISJSX1KFJ1PB34IFcNxxsOOO\nMQf3ootiWTORRiaXxMNid78u7z0RkSq99BJcfnmMqmvfHu6+G446KpZ6BhWQFJFGT3GEiEh9ygz+\nZs6EM86IgPGhh+IM1aBBChClUckl8XCtmV0EPA8sSe109zfz1isRAWDsWDjssJhi0aoV3Hcf9O9f\n6F6JiNSK4ggRkXzq0CGWNdt880hA3HJLbE+cqHm40mjkknjYEvgd0IvKIZKebItIHixfDtddFyPm\nKpK/OrNIaCvxICJNjOIIEZGGsHhxnKlavhyWLYMBA+Dpp2GHHQrdM5GcEg+HABu6+9J8d0ZE4KOP\nopbDSy/BTjtFEcnypzSYdQAAIABJREFUchWPFJEmS3GEiEhD6NEjAsalS6F16zhrtdNOcOaZcPHF\nqv0gBZVL4uE9YA3gyzz3RaTFKi2NZZi//BJuuy2KFd97L/zud/DqqyoeKSJNmuIIEZGGkFl4slu3\nSDpccUWMfLj77hhKq8BSCiCXxMMawAdm9gYrzs3UgG+RelBaCr16xeg4iMT0o49C586xreKRItLE\nKY4QEWkomYHj7bfDwQfDscfG/qIicI+REZMmKciUBpNL4uGivPdCpIWqqICRIyuTDq1awd57VyYd\nRESaAcURIiKF1KcPvPce7LknlJXFvqVLY+SDEg/SQFrVdIC7TwM+AFZLLrOSfSKyEubMgZ494Zln\nIuFQVBRTLHr2LHTPRETqj+IIEZFGoEOHqFzerl1sL18Os2bBkiXZ7ydST2pMPJjZAOB1ojjUAOA1\ns/9r787D5CqrBIy/pxuSoCKMiMtANCAwigECBjAsGuKg7KCiEJVNIEQBYRAQ1BEkoBNZBUEEWQRR\nQNawDKCQQJBmCZtsMhMlLsAAsgmYhSRn/vhu00Xb6a5AqquX9/c8/aTurVv3nrpJpU+d+33nxo6N\nDkwaqBYuhFNPhbXWgvvuK9Ptpk+HSZMc8SZp4DGPkKQ+YsyY0lTsO9+BbbaB88+H0aPh7rubHZkG\ngXqmWnwLWD8znwaIiBWB3wCXNDIwaaBpa4PLLy/FhXvuKaPezjwThg8vz2+0UXPjk6QGMY+QpL6i\ntgfEtdfC3nuX221+85tlKsZvf2vjSTVEPYWHlvZkofIsdYyUkNThttvK/+GvvlqWDzus9HaIaGpY\nktQbzCMkqS/aaqvS++HAA8vQ26OPLsnp0KEOw9USV0/h4bqIuB74ZbW8E/DfjQtJGlj+/GfYddeO\nokNrK7z97RYdJA0a5hGS1Ff9y7+Ue7hnlqkXmaXruYUHLWH1NJc8BPgJsHb1c0ZmHtrowKT+LhPO\nPrv0cnj8cVh66VJ0GDKkjH6QpMHAPEKS+oGvfAWWWaZcGWsvQjz4YLOj0gDS44iHiFgFuDYzL6uW\nl4mIEZk5q9HBSf3VE0+UKXPXXgsf/3gpQDz1VLlrkdPmJA0m5hGS1A+MGVNGOUybVpZPPBE+8hH4\n7nfh4INhqXoGykuLVs+/oF8BtW3vFlTr1m9IRFI/1dZWGgXPm1fuVjRnDpx0Euy/f7ld5qqrWnCQ\nNCiZR0hSf1DbeHKvvcooiMMPhyuuKH0gHnvMK2h6w+opPCyVmfPaFzJzXkQMaWBMUr/T1gbjxpVi\nA8DIkXDppbDGGs2NS5L6APMISepvVlwRfvUruOgimDABxo8v0zCGDbP/g96QerpKPxMR27UvRMT2\nwN8aF5LU//zoRx1FhwjYeWeLDpJUMY+QpP6oPandb7+ynAmzZ5era9JiqqfwMBH4ZkT8OSL+DHwD\nmNDYsKT+4bnn4AtfgF/8ovzf3NpaCsHjxjU7MknqM8wjJKk/23bb0niypfrqeOqpcPrppRAh1anH\nqRaZ+QfgoxHxtmr55YZHJfUD11xTpr/97W9w1FFlytuttzr1TZJqmUdIUj9X23jyQx+C004r/R8u\nuwzOOguGD292hOoH6m5PaqIglV4O110H99wDV19dbpV57bWw7rrl+U03bW58ktRXmUdIUj9W23hy\n++3hjDPg618vjc1++MMyx/jmm70Cp0XyvihSndraYLPNYO7csrzrruX/3KFDmxuXJEmS1GsiYJ99\nYPPNYY89yk9LS1k/ZIjNJ9Wleno8SIPeK6/AQQd1FB1aW+GDH7ToIEmSpEFq1VXLveS33hoWLoQF\nC0qyPG1asyNTH9Rj4SEi3hIR/xkRZ1bLq0fENo0PTeobbrsNRo2C22+HpZYqRYchQ8pIMklS98wj\nJGkAa2mBb32rdFeHUoC46SZ49tnmxqU+p54RD+cAc4H28TKPA0c3LCKpj5gzB77xjdK3Yf78UtC9\n5RaYNMkRZJK0GMwjJGkgGzOmFBuOPhomTCi9HkaOLJ3YpUo9PR4+kJk7RcR4gMz8R0REg+OSmqat\nDS64oDSNfOwx2HtvOP54WHbZ8rwFB0laLOYRkjTQ1Taf/MpXSjO0bbaBL38Zxo+Hu+6y8eQgV0/h\nYV5ELAMkQER8gHLlQhpwpk+HcePKCAcoBYeDDmpuTJLUz5lHSNJgMmpUKTR897vwX/8F55xTGk8O\nHeqw4UGsnqkWRwDXAcMj4gLgRuDQhkYlNcEjj5SCbHvRobW1o5mkJOkNM4+QpMFm6FD43vfK3S8y\nS++HOXPghhuaHZmapNvCQzUU8vfAZ4DdgV8CozNzWsMjk3rJwoVw4omw3nrw0kulcaQNJCXpzTOP\nkKRBbtddYZllyoiHTDj7bLjjjmZHpSbodqpFZmZEXJuZawF2B9GAM2sW7L576YGzzTZw5pmlr8O0\naU5Dk6Q3yzxCkga5MWPK9Ipp0+Btb4Njj4WNNoLDD4fvfKdc6dOgUE+Ph3siYv3MvKvh0Ui9oK2t\n3KHi5ZfhlFNKAfass2CPPcrj97zHgoMkLUHmEZI0mNU2ntx1VzjgADjmmHLXi69/Hf7yF6/4DQL1\nFB42BL4YEX8CXgGCchFj7YZGJjVAW1tpHjlnTlled1247DIYMaKpYUnSQGYeIUkqllsOzj0Xdtih\n3PFil13Klb9hw2w8OcDV01zyU8AHgHHAtsA21Z89iogtIuLRiJgZEYd18fzQiLioev6OiBhR89zh\n1fpHI+JTPe0zIlap9jGz2ueQTsf6bERkRIyulkdExOyIuK/6Ob2e96T+7ZRTOooOLS2w444WHSSp\nwcwjJEmvt8MO8NWvlseZMHs2XHppc2NSQ/VYeMjMPwHLU5KEbYHlq3XdiohW4FRgS2BNYHxErNlp\nsz2B5zNzNeBEYHL12jWBnYEPA1sAp0VEaw/7nAycWO3r+Wrf7bEsCxwAdO5k8ofMHFX9TOzpPan/\nev55+OIX4Ze/LEXV1tbSbHezzZodmSQNbOYRkqQubb11aTzZUn0lPfVU+MlPSiFCA06PhYeIOAC4\nAHhX9fPziNi/jn1vAMzMzD9m5jzgQmD7TttsD/ysenwJ8ImqA/b2wIWZOTczHwNmVvvrcp/Va8ZV\n+6Da5w41x5lESSjm1BG3Bpjrr4eRI+Hii+Goo0ojyUmTHM0lSb3BPEKS1KX2xpNHHw2XXw6bbAIT\nJ8JWW8Hjjzc7Oi1h9fR42BPYMDNfAYiIyUAbcEoPr1sJ+EvN8l8p8zy73CYz50fEi8AK1frbO712\npepxV/tcAXghM+d33j4i1gOGZ+Y1EXFIp+OvEhH3An8Hvp2Z03t4T+on2trKbYLvv7/8P7bmmjBl\nCnzkI+X5TTdtbnySNIiYR0iSulbbeHK77eDHP4ZDDoG11oLTToP3v9/bzQ0Q9RQeAlhQs7ygWtfn\nRUQLcALl3uGdPQm8LzOfjYiPAFdExIcz8++d9jEBmADwvve9r8ERa0loaytTKObOLcvjx5dbBg8b\n1ty4JGmQMo8wj5CknrW0wL77wuabl7tfjB9f5kdDue2mw5X7tXqaS54D3BERR0bEkZQrCGfV8brH\ngeE1yytX67rcJiKWApYDnu3mtYta/yywfLWP2vXLAiOBaRExC/goMCUiRlfDL58FyMy7gT8Aa3R+\nE5l5RmaOzszRK664Yh1vW800d265LXB70aGlpRRMLTpIUtOYR5hHSFL91lgDbr0VPvlJWLCg/Myd\nW0Y+qN+qp7nkCcAewHPVzx6ZeVId+74LWL3qEj2E0uRpSqdtpgC7VY93BG7KzKzW71x1q14FWB24\nc1H7rF4ztdoH1T6vzMwXM/OdmTkiM0dQkp3tMnNGRKxYNZkiIlatjvHHOt6X+qjf/Q422KD0cGht\n7WggOXZssyOTpMHLPEKStNiWWgqOPLIk8wALF8Idd8DLLzc1LL1xi5xqERHDgInAasADwGk1cx97\nVM213A+4HmgFzs7MhyLiKGBGZk6hXPE4PyJmUpKRnavXPhQRFwMPA/OBfTNzQRXXP+2zOuQ3gAsj\n4mjgXnq+mvIx4KiIeBVYCEzMzOfqfX/qOxYsgOOOg//8T3jHO+Cqq2CFFZwOJknNZB4hSXpTxoyB\nqVPhN7+B3/++3J5unXXgvPNg442bHZ0WU+QiblcSERcBrwLTKbedmpWZB/ZibH3O6NGjc8aMGc0O\nQ5W2tnK73xtugAcegM9+Fk4/Hd75zmZHJkn9U0TcnZmjl9C+zCM6MY+QpDdh+nTYbTeYNQsOPRS2\n3BJuu80rjX1Id3lEd80l18zMtaodnEUZoij1Ce3/x7z6alk+4ojyE/2iXZkkDQrmEZKkJWfTTcst\n6w4+GCZPhh/8oCT/Q4faeLIf6K7Hw6vtDxZnaKTUaE8+CXvs0VF0aO/lYNFBkvoU8whJ0pK17LLw\nk5+UkQ+ZpffDnDlw003Njkw96K7wsE5E/L36eQlYu/1xRPy9m9dJDXPxxTByZBlhtfTSpegwZIgN\nJCWpDzKPkCQ1xj77dNy2LhMuughmzmxuTOrWIqdaZGZrbwYidee552C//UpPmQ02KD1lnnvOBpKS\n1FeZR0iSGmbMmDLKYerUMgz6pJNK48njjy9FCYdC9znd9XiQmq6tDc44A66+Gl54AY46Cg4/vNxh\nByw4SJIkSYPSmDEdXwb23BO+/GX4ylfgyivLnw895BXKPsTCg/qsm26CT30K5s8vRcuzziq9HSRJ\nkiTpNSuvDNddBz/+MRx0UHnc0mLjyT6kux4PUtO0tcHnPleKDlD+3/i//2tuTJIkSZL6qJYW2Hff\nMj8bOhpPXnttc+MSYOFBfcy8efCtb8Emm5TmkUOH2kBSkiRJUp123BGWWaYMmc4soyCuu67ZUQ16\nFh7UZzz4IGy4IXzve+UOOf/zP6VfzKRJjpCSJEmSVIcxY8qXh2OOgXPOgXe/G7bcsvR9eOWVZkc3\naNnjQU3V1lZ6OTz5JJx5Jiy/fOkHs9125fnanjGSJEmS1KPaLxE77wzf/jaccAL8+tdw6KHw7LM2\nnuxlFh7UNG1tMG5cmXoF8LGPwSWXwIorNjcuSZIkSQPEsGFw3HGw7balCNF+u81hwxxW3YucaqGm\nyIRjj+0oOrS0lDtYWHSQJEmStMR9/OMwYUJH74fZs+HCC5sd1aBh4UG97sknS8Hx8stLwaG1tTSR\n3GyzZkcmSZIkacDaYosy0qGl+hp82mkweTIsWNDcuAYBp1qoV/3qVzBxIvzjH3DyybDeenDLLU6x\nkiRJktRg7Y0np02DddaBs86Cww6Dq66Cn/0MPvCBZkc4YFl4UMO1tcF//zfcfnvp57L++nDeefDB\nD5bnN964ufFJkiRJGiRqG09uuSVccAHst18pRBx/PKy1Ftx8s1dGlzALD2qotrYyhWLu3LK8117l\nVrpL+S9PkiRJUjNFwJe+VPo/fPnLZWh2S0tZP2SIzSeXIHs8qGFefhkOPLCj6NDaCquuatFBkiRJ\nUh8yfDhcf31pRLdwYen5MHcuTJ3a7MgGDAsPaohbb4VRo+DOO0uhobW1FA3Hjm12ZJIkSZLUSUsL\nHH54aT4JpQBx443wt781N64BwmvPWqLmzIHvfKfcKvf97y/To5ZeuvRvcZqUJEmSpD5rzBi46aby\n89e/luaTI0fCT38K22zT7Oj6NQsPWiLa2uAXv4BrroHHHoO99y69WZZdtjxvwUGSJElSn1fbfHLi\nRNhllzIFY889Yeed4a67vKL6Blh40Js2fTqMGwfz55fl446Dr3+9uTFJkiRJ0puyzjql0HDkkTB5\nMpx9dmk8OXSojScXkz0e9KY88giMH99RdGhthXnzmhuTJEmSJC0RQ4fC978PEyZAZun9MGcO3HBD\nsyPrVyw86A1ZsABOOAHWXRdeeqk0jrSBpCRJkqQBabfdYJllyoiHzNL/4fbbmx1Vv2HhQYvtD3+A\nzTYr0yk+9Sl49NHSPHLSJEccSZIkSRqAxowpX3aOOQZOPrkUIDbeGL75zXLrTXXLHg+qS1tbuY3t\nCy/AaaeV0Q3nngu77lo+c+95jwUHSZIkSQNYbePJ3XaDgw4q0zCuvrpclX3iCRtPLoKFB/Wora00\nj5wzpyyvvz5ceikMH97cuCRJkiSpKd7+9nKbzU9/ulyN3X33ckV22DCHgXfBqRbqVma5LWZ70aGl\nBXbYwaKDJEmSJLH11vDVr5bHmTB7Nlx0UXNj6oMsPGiRnngCttuujG5oaSnTK4YOLf0dJEmSJEnA\nVluVxpMt1dfrU08tt99csKC5cfUhTrXQP8mEX/wC9t+/jHQ46SQYPRpuucUpS5IkSZL0Ou2NJ6dN\ng7XXhrPPhsMOgyuuKI3x/u3fmh1h01l40Gva2kpflFtvLUWGMWPK52SNNcrzG2/c1PAkSZIkqW+q\nbTy51VZw4YWw774walS5E8YGG8D06YP2Sq6FBwFw221lCsW8eWV5333hhz8s0yskSZIkSXWKgPHj\nS5Fhn33KHS/ap2EMHToom0/a40E89RRMmNBRdGhthZVWsuggSZIkSW/Ye98LV14JO+4ICxeWn7lz\nYerUZkfW6yw8DGKZZQTQhz8Mjz4KSy9dig1DhpTinCRJkiTpTYiAgw4qt9mEUny46CL43/9tbly9\nzKkWg9RTT5W7vlx2WZludM458OKLpR/KIJ12JEmSJElL3pgxcNNNZaTDnDlwyimlCeX3vgdf+9qg\nGGpu4WGQue02OPlkuPbaMrVi8uRSgFuq+pdgwUGSJEmSlrDa5pMTJ5beDwcdBJdcUm4n+NhjA/oK\nsIWHQeSqq2CHHcronpYW+PnPS88TSZIkSVIv+dd/hSlTyheyr361fCmLKNMxBmjjSXs8DAKZ5baY\nn/98KTpA+Xc9a1Yzo5IkSZKkQSoCdtmljHaA8qVt9uzS/2EAsvAwwP3pT7DllrDHHrDGGqWIZgNJ\nSZIkSeoDtt0Wllmm43abP/oRHHVUxy0HBwgLDwPUwoVw6qkwciTcemvpX3LvvaWnyaRJA3YEjyRJ\nkiT1H2PGlC9nRx8N11wDn/scHHEErL8+3H13s6NbYuzxMMC0tZX+JDfeCPffD5tvDmecASNGlOdr\ne5pIkiRJkpqs9kvaVluVng8TJ8KGG5bHq69evtj14y9yDR3xEBFbRMSjETEzIg7r4vmhEXFR9fwd\nETGi5rnDq/WPRsSnetpnRKxS7WNmtc8hnY712YjIiBjd0zH6q+nT4eMfhxNOKEWHb34Trr++o+gg\nSVJ/Yh4hSRqUttsOHn64zJn/+c/LCIjNNitXmfuphhUeIqIVOBXYElgTGB8Ra3babE/g+cxcDTgR\nmFy9dk1gZ+DDwBbAaRHR2sM+JwMnVvt6vtp3eyzLAgcAd9Ss6/IYS+4M9K6774addoJXXy3Lra3w\ntreVniWSJPU35hGSpEFt+eVho406ej/MnQsHHAAvvNDcuN6gRo542ACYmZl/zMx5wIXA9p222R74\nWfX4EuATERHV+gszc25mPgbMrPbX5T6r14yr9kG1zx1qjjOJklDM6XTsro7Rr7zyChx8MGywQek/\nMmSIzSMlSQOCeYQkaXAbOxaGDi1f8JZaCmbMgDXXhMsua3Zki62RhYeVgL/ULP+1WtflNpk5H3gR\nWKGb1y5q/QrAC9U+XnesiFgPGJ6Z17yB+IiICRExIyJmPPPMM929317T1gbf/z6cdBKstRYcfzzs\nvTfMnAnTptk8UpI0IJhHSJIGt/bGk5MmwS23wF13wbvfDZ/9LHz60/D44x1fDvv4NIwB3VwyIlqA\nE4Dd3+g+MvMM4AyA0aNH55KJ7I1ra4Nx42BOdc1l+HC4+Wb42MfKss0jJUlaMgZiHiFJ6mc6f8G7\n80448cTS92GNNcpc+4ULy5D3Pnz1uZEjHh4Hhtcsr1yt63KbiFgKWA54tpvXLmr9s8Dy1T5q1y8L\njASmRcQs4KPAlKoxVD3x9SmZZXRDe9Ehoox0aC86SJI0gJhHSJLU2dJLw6GHwgMPwHveUwoPCxaU\neffTpjU7ukVqZOHhLmD1qkv0EEoDpimdtpkC7FY93hG4KTOzWr9z1a16FWB14M5F7bN6zdRqH1T7\nvDIzX8zMd2bmiMwcAdwObJeZM7o5Rp/0yCNlis+ll5b+Iq2tMGwY/Pu/NzsySZIawjxCkqRFWW01\nOP/8MtIBSvHhoYdKE8A+qGFTLTJzfkTsB1wPtAJnZ+ZDEXEUMCMzpwBnAedHxEzgOUoCQLXdxcDD\nwHxg38xcANDVPqtDfgO4MCKOBu6t9t1dfIs8Rl8yezYccwz84AflLhVnngkf+lCZ4jN2bJ8dSSNJ\n0ptiHiFJUg822qiMcrj2WrjvPrjgArj1VvjRj2CbbZod3etEKfKrHqNHj84ZM2b0yrHa2kqR4YYb\nSs+QXXaB446Dd72rVw4vSRqEIuLuzBzd7DgGqt7MIyRJg9D06TBxIjz8MHzmM7DrruVxL12x7i6P\nGNDNJfurKVPKv5MFC0ofh5NPhv33b3ZUkiRJkqQ+a9NN4d57S2PAI48st92MKHP0m9x4spE9HrSY\n5s2DY4+Fz32uFB2g9HN4+eXmxiVJkiRJ6geGDIHDD4evfa0sZ5b5++ec09SwLDz0ETfeCOusUxqU\nrr9+KUq1tpZ/N2PHNjs6SZIkSVK/8ZnPwDLLlCvZEWUe/xe+AE880ZRwnGrRRG1tcMUVMGMG3HQT\nrLoqXH01bL11eW7aNBtISpIkSZIW05gx5er2tGnl8dSpMHkyXHUVHHEEbLhhaUTZS184LTw0yc03\nw+abl9uuAuy1F5xyShnpAOXv3oKDJEmSJOkNqf1SOXZsaTZ54IFwyCFlFEQEDB3aK/0fnGrRJOed\n11F0aG0tox3aiw6SJEmSJC1RH/hAGfGw666l98PChaXR4LRpDT+0hYcm2WuvUlyyj4MkSZIkqddM\nnFj6P/Til1GnWjRJ+zQb+zhIkiRJknpNbf8HezwMfPZxkCRJkiT1ul7+MupUC0mSJEmS1DAWHiRJ\nkiRJUsNYeJAkSZIkSQ1j4UGSJEmSJDWMhQdJkiRJktQwFh4kSZIkSVLDWHiQJEmSJEkNY+FBkiRJ\nkiQ1jIUHSZIkSZLUMJGZzY6h34iIZ4A/LeHdvhP42xLe50DjOeqZ56hnnqOeeY56NtDP0fszc8Vm\nBzFQNSiP6GsG+mdkSfAc9cxz1DPPUc88Rz1b0udokXmEhYcmi4gZmTm62XH0ZZ6jnnmOeuY56pnn\nqGeeI6l7fkZ65jnqmeeoZ56jnnmOetab58ipFpIkSZIkqWEsPEiSJEmSpIax8NB8ZzQ7gH7Ac9Qz\nz1HPPEc98xz1zHMkdc/PSM88Rz3zHPXMc9Qzz1HPeu0c2eNBkiRJkiQ1jCMeJEmSJElSw1h46CUR\ncXZEPB0RD9ase0dE/Doi/rf681+aGWOzRcTwiJgaEQ9HxEMRcUC13vNUiYhhEXFnRNxfnaPvVutX\niYg7ImJmRFwUEUOaHWuzRURrRNwbEVdXy56jGhExKyIeiIj7ImJGtc7PWo2IWD4iLomI30fEIxEx\nxnMkdTC36Z55TX3MbepjXtMzc5ueNTO3sfDQe84Ftui07jDgxsxcHbixWh7M5gNfz8w1gY8C+0bE\nmnieas0FxmXmOsAoYIuI+CgwGTgxM1cDngf2bGKMfcUBwCM1y56jf7ZZZo6quY2Sn7XX+yFwXWZ+\nEFiH8u/JcyR1OBdzm+6Y19TH3KY+5jX1MbfpXtNyGwsPvSQzbwGe67R6e+Bn1eOfATv0alB9TGY+\nmZn3VI9fonwQVsLz9JosXq4Wl65+EhgHXFKtH9TnCCAiVga2Bn5aLQeeo3r4WatExHLAx4CzADJz\nXma+gOdIeo25TffMa+pjbtMz85o3xc9bpdm5jYWH5np3Zj5ZPf4/4N3NDKYviYgRwLrAHXieXqca\nancf8DTwa+APwAuZOb/a5K+UxGYwOwk4FFhYLa+A56izBG6IiLsjYkK1zs9ah1WAZ4BzqqGtP42I\nt+I5knriZ6QL5jXdM7fpkXlNfcxtutfU3MbCQx+R5fYi3mIEiIi3AZcCB2bm32uf8zxBZi7IzFHA\nysAGwAebHFKfEhHbAE9n5t3NjqWP2yQz1wO2pAz//Vjtk37WWApYD/hxZq4LvEKnoYeeI6l7fkYK\n85qemdssmnnNYjG36V5TcxsLD831VES8F6D68+kmx9N0EbE05ZfzBZl5WbXa89SFamjUVGAMsHxE\nLFU9tTLweNMCa76Nge0iYhZwIWUo4g/xHL1OZj5e/fk0cDkl0fOz1uGvwF8z845q+RLKL2vPkdQ9\nPyM1zGsWj7lNl8xr6mRu06Om5jYWHpprCrBb9Xg34MomxtJ01Xy1s4BHMvOEmqc8T5WIWDEilq8e\nLwNsTpkzOhXYsdpsUJ+jzDw8M1fOzBHAzsBNmflFPEeviYi3RsSy7Y+BTwIP4mftNZn5f8BfIuLf\nqlWfAB7GcyT1xM9IxbymPuY23TOvqY+5Tc+andtEGU2hRouIXwJjgXcCTwFHAFcAFwPvA/4EfD4z\nOzdpGjQiYhNgOvAAHXPYvkmZD+l5AiJibUrTl1ZK4fDizDwqIlalVMHfAdwLfCkz5zYv0r4hIsYC\nB2fmNp6jDtW5uLxaXAr4RWYeExEr4GftNRExitLIawjwR2APqs8dniPJ3KYH5jX1Mbepn3nNopnb\n1KeZuY2FB0mSJEmS1DBOtZAkSZIkSQ1j4UGSJEmSJDWMhQdJkiRJktQwFh4kSZIkSVLDWHiQJEmS\nJEkNY+FB6oMiYkFE3BcRD0XE/RHx9YhoqZ4bHREnV4+HRsRvqm13iohNq9fcV90Lu1HxjY2IjWqW\nJ0bErj285siIOLhRMTVLRFxS3cKJiJgVEdM7PX9fRDxYPR4bES9W6x6JiCN62PeKEXFd46KXJKlr\nNbnIgxHxq4h4S7Nj6kuiuCki3l4tZ0T8vOb5pSLimYi4ulrevVq+LyIejoi9e9j/WhFxbkPfhNSL\nLDxIfdPszBw+fmLvAAAIXUlEQVSVmR8GNge2pNwfncyckZlfq7Zbt1o3KjMvAr4IfL9ant3TQapf\nmm/k/4GxwGuFh8w8PTPPewP7qUtELNWofb+ZY0XEh4HWzPxjzeplI2J49fyHunjZ9MwcBYwGvhQR\n6y1q/5n5DPBkRGxcb0ySJC0h7bnISGAeMLGeF72J3GJR++uTOQCwFXB/Zv69Wn4FGFlz4Wdz4PFO\nr7moygHGAt+LiHcvaueZ+QCwckS8bzFikvosCw9SH5eZTwMTgP2qX+ZjI+LqiHgX8HNg/ap6vg/w\neWBSRFwAEBGHRMRdEfG7iPhutW5ERDwaEecBDwLDI+KTEdEWEfdUVzXeVm07KyK+W61/ICI+GBEj\nKMnHf1TH3bR2NENE7F0d8/6IuLSnKyQRcW5EnB4RMyLifyJim2r97hExJSJuAm7s5v28NSKuqY73\nYETsVK3/r+qKwu8i4riaY+1Yc+yXqz/HRsT0iJgCPFyt+1JE3Fm9x59ERGsX4X8RuLLTuouBnarH\n44FfLuLv9RXgbmC16u9kenWe74ma0STAFdVxJElqlunAagARcVD1+/bBiDiwWld3blErIqZFxA+j\nY2TFBtX6IyPi/Ij4LXB+RLRGxLE1OcA+1XbvjYhbal6/abXtudXyAxHxHzXHGl09fmdEzKoe15Vv\ndKGrHOBaYOvqcXc5wNPAH4D3R8QG1Xm6NyJui4h/q9n0KmDnRRxf6lcsPEj9QHVFvRV4V826p4G9\nqK6gZ+ZPgCnAIZn5xYj4JLA6sAEwCvhIRHysevnqwGnViIpXgG8D/56Z6wEzgINqDv+3av2PgYMz\ncxZwOnBiddzXTS0ALsvM9TNzHeARYM863uKIKs6tgdMjYli1fj1gx8z8eDfvZwvgicxcp7oqc11E\nrAB8GvhwZq4NHF1HDOsBB2TmGlFGKuwEbFxdmVhA11/+N6YUD2pdCnymerwtJWn4J1WMHwUeAp4G\nNq/O807AyTWbzgA2rSN+SZKWuCijALYEHoiIjwB7ABtSfoftHRHrVpsuTm5R6y3V79qvAmfXrF+z\nev14Si7xYmauD6xfHXcV4AvA9dXr1wHuo+QIK2XmyMxcCzinjrdZT77RWVc5wIXAzlUeszZwR1cH\nizJFc1VgJvB7YNPMXBf4DvC9mk3NATRg9NrQJUm97pPVz73V8tsov0j/DPwpM2+v1n+U8sv9txEB\nMARoq9nPZdWfd9Pxhbo7IyPiaGD56pjX1/GaizNzIfC/EfFH4IPV+l9n5nM9vJ/pwPERMRm4OjOn\nV0nSHOCsKHMrr64jhjsz87Hq8SeAjwB3VedkGUpxoLP3As90Wvcs8HxE7EwpvPyj0/ObRsS9wELg\nvzLzoYhYDvhRRLQXOdao2f5p4F/riF+SpCVpmYi4r3o8HTgL+ApweTVqj4i4jPLFeAqLl1vU+iVA\nZt4SEW+PiOWr9VNqpo1+Eli7ZtTicpQc4C7g7IhYGrgiM++r8ohVI+IU4Brghjreaz35xi2dXvOO\nzHypdkVm/i7KyNDxlNEPne0UEZsAc4F9MvO5KNMzfxYRqwMJLF2zvTmABgwLD1I/UFXGF1B+AXXV\nN6DLl1H6Pfyk075GUK5E1G736+qKQlfmVn8uoL7/M84FdsjM+yNid8o8xp7kIpY7x/lP7wcgSp+E\nrYCjI+LGzDyqGq75CWBHYD9gHDCfaqRXlPmnQ2p20/lYP8vMw3uIezYwrIv1FwGnArt38dz0zNym\n07r/AJ6iXK1poRRN2g2rjiNJUm+aXY0keE1VRFiUxcktatWbA+yfmf90MaMajbA1cG5EnJCZ50XE\nOsCnKFNDPw98mZocgH/+3V1XvtHJ/IhoqS6c1JoCHEfJf1bo9NxFmblfp3WTgKmZ+ekqR5tW85w5\ngAYMp1pIfVxErEiZ2vCjzOz8y7k71wNfjo5+DStF6QvR2e3AxhHRPnfzrRGxRhfb1XoJWHYRzy1L\naYi4NPX3JvhcRLRExAcoQw8f7WKbLt9PRPwr8I/M/DlwLLBetc1ymXkt5Uv9OtU+ZlFGMgBsx+uv\nKtS6Edix/XxFxDsi4v1dbPcI1ZzXTi4HfkB9oz2gXLl5skpedqFMq2m3BmW+rCRJzTYd2CEi3hIR\nb6VMa+w85RIWL7do7820CWU6xYtdbHM98JUqtyAi1qj2+X7gqcw8E/gpJQd4J9CSmZdSpnu0N3Ge\nRUcOsCOLVm/+9CglZ+nsbOC7VXPIeixHRxPK3Ts9Zw6gAcMRD1Lf1D68cWlKhf584ITF2UFm3lD1\nKmirrlC8DHyJMnKhdrtnqpEJv4yIodXqbwP/083urwIuiYjtgf07PfeflDmNz1R/LqpAUevPwJ3A\n24GJmTmn81WVbt7PasCxEbEQeJUyDHRZ4MpqjmXQMa/0zGr9/cB1vP4KR+2xHo6IbwM3VCMjXgX2\nBf7UadNrKFc0ftPp9S8Bk6HHq0PtTgMujXJL0s5xbVYdR5KkpsrMe6Lc4vHOatVPM/Pe6kp97XaL\nk1vMqaYgLk0ZmdCVn1L6Qd0T5RfrM8AOlN/Bh0TEq5S8YFdgJeCc6LizRvvoxeOAiyNiAt38Xu0m\n3+g85bI9B5jZ6fV/5fW9mnryA8pUi293EZc5gAaMWLwLqJK0ZFUJzNWZeUmzY1lcUW6ZNZXShHJB\nT9u/wWPcAmyfmc83Yv+SJDVLREyjNK6e0exYFldEvBc4LzM3b9D+hwI3A5tk5vxGHEPqTU61kKQ3\nqGp6dQTl6soSV02zOcGigyRJfUtmPgmcGRFvb9Ah3gccZtFBA4UjHiRJkiRJUsM44kGSJEmSJDWM\nhQdJkiRJktQwFh4kSZIkSVLDWHiQJEmSJEkNY+FBkiRJkiQ1jIUHSZIkSZLUMP8PAwrFicfx4ksA\nAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 1080x360 with 2 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "h_PG2SNU-lU1",
        "colab_type": "text"
      },
      "source": [
        "**Conclusion**: In production, while pore pressure decreases, differential pressure increases, the **formation compressibility increases**, due to pore collapse. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WNtaahQ-4W6K",
        "colab_type": "text"
      },
      "source": [
        "# ***"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Xkuzd15qkyy8",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        },
        "outputId": "2c94902c-248b-4d5e-f3e9-be719a477e7d"
      },
      "source": [
        "# psi to GPa\n",
        "C_formation = 5E-6 / 6.89E-6 # psi^-1 to GPa^-1\n",
        "C_rock = C_formation / ((1 - poro) / poro)\n",
        "K_rock = 1 / C_rock\n",
        "K_rock"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "1.378"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 46
        }
      ]
    }
  ]
}