files/assignment1_solution-2018.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem 2 (30 points)\n",
    "\n",
    "Determine the principle stresses and directions for the following stress tensors:\n",
    "\n",
    " a. $$\\mathbf{S} = \\begin{bmatrix} 18 & 0 & 24 \\\\ 0 & -50 & 0 \\\\ 24 & 0 & 32 \\end{bmatrix}$$\n",
    " \n",
    " b. $$\\mathbf{S} = \\begin{bmatrix} 3 & -10 & 0 \\\\ -10 & 0 & 30 \\\\ 0 & 30 & -27 \\end{bmatrix}$$\n",
    " \n",
    " c. $$\\mathbf{S} = \\begin{bmatrix} 2 & -1 & 1 \\\\ -1 & 0 & 1 \\\\ 1 & 1 & -2 \\end{bmatrix}$$\n",
    " \n",
    " Work (a) by hand for practice.  You can use a computer to solve (b) and (c) if you wish.\n",
    " \n",
    " **Solution:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eigenvalues: \n",
      " [  0.  50. -50.]\n",
      "Eigenvectors: \n",
      " [[-0.8 -0.6  0. ]\n",
      " [ 0.   0.   1. ]\n",
      " [ 0.6 -0.8  0. ]]\n"
     ]
    }
   ],
   "source": [
    "import scipy.linalg\n",
    "import numpy as np\n",
    "\n",
    "values, vectors = scipy.linalg.eig(np.array([[18,0,24],[0, -50, 0],[24, 0, 32]]))\n",
    "print(\"Eigenvalues: \\n\", np.real(values))\n",
    "print(\"Eigenvectors: \\n\", vectors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eigenvalues: \n",
      " [  2.30000000e+01   2.66453526e-15  -4.70000000e+01]\n",
      "Eigenvectors: \n",
      " [[ 0.3940552   0.9124486  -0.11026357]\n",
      " [-0.78811041  0.27373458 -0.55131785]\n",
      " [-0.47286624  0.30414953  0.82697677]]\n"
     ]
    }
   ],
   "source": [
    "values, vectors = scipy.linalg.eig(np.array([[3,-10,0],[-10, 0, 30],[0, 30, -27]]))\n",
    "print(\"Eigenvalues: \\n\", np.real(values))\n",
    "print(\"Eigenvectors: \\n\", vectors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eigenvalues: \n",
      " [-2.77845712  2.48928857  0.28916855]\n",
      "Eigenvectors: \n",
      " [[ 0.26826216  0.93698901  0.22380127]\n",
      " [ 0.41025777 -0.32131516  0.85348997]\n",
      " [-0.87162146  0.13714285  0.47060372]]\n"
     ]
    }
   ],
   "source": [
    "values, vectors = scipy.linalg.eig(np.array([[2,-1,1],[-1, 0, 1],[1, 1, -2]]))\n",
    "print(\"Eigenvalues: \\n\", np.real(values))\n",
    "print(\"Eigenvectors: \\n\", vectors)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem 3 (40 points)\n",
    "\n",
    "The following table contains the estimated bulk mass densities as a function of depth for an offshore location in Brazil. Water depth is 500m. Measurements indicate that porosity of shale layers estimated through resistivity measurements.\n",
    "\n",
    "|Depth (m)| Bulk mass density (kg/m$^3$) | Shale porosity |\n",
    "|:-------:|:-------------------------:|:--------------:|\n",
    "| 0   | 1025 | | \n",
    "| 100 | 1026 | |\n",
    "| 200 | 1026 | |\n",
    "| 300 | 1030 | |\n",
    "| 400 | 1030 | |\n",
    "| 500 | 1031 | |\n",
    "| 600 | 1900 | |\n",
    "| 700 | 2190 | |\n",
    "| 800 | 2200 | |\n",
    "| 900 | 2230 | |\n",
    "| 1000| 2235 | |\n",
    "| 1100| 2240 | |\n",
    "| 1200| 2275 | 0.305 |\n",
    "| 1300| 2305 | 0.297 |\n",
    "| 1400| 2310 | 0.286 |\n",
    "| 1500| 2308 | 0.281 |\n",
    "| 1600| 2310 | 0.285 |\n",
    "| 1700| 2305 | 0.293 |\n",
    "| 1800| 2310 | 0.307 |\n",
    "| 1900| 2324 | 0.305 |\n",
    "| 2000| 2319 | 0.298 |\n",
    "\n",
    "\n",
    " a. Plot $S_v$ as a function of depth (in SI units)\n",
    "\n",
    " b. Plot hydrostatic water pressure as a function of depth. Assume the density of brine water is 1031 kg/m$^3$\n",
    "in the rock pore space (in SI units).\n",
    "\n",
    " c. Additional compaction lab measurements on shale cores indicate a good fit of the porosity-effective vertical stress relation through the equation $\\phi = \\phi_0 \\exp(-\\beta (S_v - P_p))$, with parameters $\\phi_0=0.38$ and $\\beta = 3 \\cdot 10^{-2}$ MPa$^{-1}$. Estimate the pore pressure in the shale. Is there overpressure? At what depth does it start?\n",
    "\n",
    " d. Plot vertical effective strees, $\\sigma_v^{eff}$ as a function of depth (in SI units)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3Xl8DPcfx/HXF3GGUDTOuinqqKBF\nD3FrHRVn1E1D66xWq3W26qoeKG3dRSVRVx0l6khc7Q+J+xbaEkeJOhJy5/v7Y5aGknt3NtnP8/HY\nh93ZnZ33zsZ8dma+8/0qrTVCCCEcVxazAwghhDCXFAIhhHBwUgiEEMLBSSEQQggHJ4VACCEcnBQC\nIYRwcFIIhBDCwUkhEEIIByeFQAghHFw2swMkR6FChXTp0qVTNe/du3fJkydP+gZKB5IrZSRXykiu\nlLHXXJC2bEFBQaFa68JJvlBrbfc3Nzc3nVr+/v6pnteaJFfKSK6UkVwpY6+5tE5bNiBQJ2MbK4eG\nhBDCwUkhEEIIByeFQAghHFyGOFn8ODExMYSEhBAZGZno61xcXDh58qSNUiVfWnLlzJmTEiVK4OTk\nlM6phBCOKMMWgpCQEPLmzUvp0qVRSj3xdWFhYeTNm9eGyZIntbm01ty4cYOQkBDKlCljhWRCCEeT\nYQ8NRUZGUrBgwUSLQGaklKJgwYJJ7gkJIURyZdhCADhcEbjPUT+3EMI6MnQhEEKIzOrOnTsMHTqU\n8PBwqy9LCoEQQtiZEydOUKdOHWbPns3Ro0etvrwMe7LYXly5coV3332X4OBg7ty5Q/HixfH39zc7\nlhAig1q5ciW9evXC2dmZ7du3Ex8fb/Vlyh5BGnXv3p127doRGBjImTNnmDlzptmRhBAZUGxsLB98\n8AEdO3akWrVqBAUF8corr9hk2VII0iAuLo6AgABeffXVB9OqVatmYiIhREZ0/fp1mjdvzrRp03jn\nnXfYsWMHxYsXt9nypRCkQdasWWnSpAk1atSgf//+7Nmzx+xIQogMJjAwEDc3N/bs2cOiRYuYPXs2\n2bNnt2mGTHGOYNiwYRw6dOixz8XFxZE1a9YUv2fNmjWZPn16kq/btGkTe/bsYd26dbRo0YKlS5fy\nxhtvpHh5QgjHs3DhQt555x2KFCnCb7/9Rq1atUzJkSkKgZmUUrz00ku89NJL3Lx5kyNHjnD58mWq\nVKlCw4YN6dOnDzNnzsTZ2dnsqEIIOxEVFcWQIUOYO3cuTZo0wcfHh0KFCpmWJ1MUgsR+uVuzi4nN\nmzfj7u5O9uzZuXbtGrt372bhwoXExsZy5MgRsmTJQt26daUICCEeCAkJoX379uzbt4+PPvqICRMm\npOqoRXrKFIXALCtXruSdd97B2dmZHDlyMGHCBOrVq8etW7fw9fUlMDCQefPmmR1TCGEnAgIC6NSp\nExEREaxatQoPDw+zIwFSCNLkSRv5/Pnzs2PHDiZMmEC2bLKKhXB0WmumT5/OiBEjqFChAmvWrOHZ\nZ581O9YDVttKKaVyAjuBHJblrNRaj1NK/QC8Cty2vLSX1vrxZ3ozsGPHjpkdQQhhB+7evUvfvn1Z\nvnw5Hh4e/PDDD3bXI7I1f65GAY201uFKKSdgt1Jqk+W5EVrrlVZcthBCmC44OJh27dpx4sQJpkyZ\nwgcffGCXnUZarRBYBk6+31uSk+WmrbU8IYSwJxs2bKBbt25ky5aNzZs306RJE7MjPZFVLyhTSmVV\nSh0CrgFbtNZ7LU9NVEodUUp9rZTKYc0MQghhS3FxcYwdO5bWrVtTvnx5goKC7LoIACjjh7uVF6JU\nfmANMBi4AVwFsgNzgXNa608fM48X4AXg6urq5uvr+9DzLi4ulC9fPsllp/aCMmtLa67g4GBu376d\n9AtTKDw83C6bu0qulJFcKZNeue7cucOkSZPYu3cvLVu2ZNiwYWm+Sjgt2dzd3YO01rWTfKHW2iY3\nYBzw/iPTGgIbkprXzc1NP+rEiRP/mfY4d+7cSdbrbC2tuZL7+VPK39/fKu+bVpIrZSRXyqRHroMH\nD+oyZcpoJycnPWfOHB0fH5/2YDpt2YBAnYzts9UODSmlClv2BFBK5QKaAKeUUkUt0xTwBiDNa4QQ\nGdrSpUupV68e0dHR7Nq1Cy8vL7s8Kfwk1mw1VBRYrJTKinEu4iet9Qal1HalVGFAAYeAAVbMIIQQ\nVhMdHc17773HrFmzaNiwIb6+vri6upodK8Ws2WroCPD8Y6Y3stYyhRDCVi5fvkzHjh357bffeO+9\n95gyZUqGvYA0Y6YWQggT7dq1i06dOhEWFoavry+dO3c2O1KayHgEaXTlyhW6dOlC7dq1qVixIu7u\n7mZHEkJYidaamTNn0qhRI/LmzcvevXszfBEA2SNIs+7du/PWW289+GOwxUDTQgjbu3fvHm+99Rbe\n3t60adOGJUuW4OLiYnasdCF7BGkgQ1UK4RjOnTtHvXr18PHx4bPPPmPNmjWZpghAJtkjGDYMnjBA\nGXFxuUjNdVs1a0JSA5QlHKryjTfeoEePHjRo0CDlCxNC2K2NGzfy5ptvopRi06ZNNG/e3OxI6U72\nCNJo06ZNrFq1ChcXF1q0aMHPP/8MQL9+/ZI1/+bNm1m6dOlD07QNrvYWQiQuPj6e8ePH8/rrr1O6\ndGmCgoIyZRGATLJHkNgv97CwCKt2+fq4oSqbNWtGcHAwo0aN4sSJE7Rp04YiRYrQokUL+vbty9df\nf82IESMoVKgQe/fuZfbs2Vy9ehUPDw/atGlDz549mTZtGlFRURQoUIDPPvvMavmFEP/1zz//0L17\ndzZu3EjPnj357rvvyJUrl9mxrEb2CNJg8+bNREdHAzwYqrJp06YcOHCAli1bMnHiRPLkycMrr7zC\n7t27WbBgAZ07d+aHH36ga9euTJo0iaioKCpVqsTBgwfp0qULI0eOZOnSpXTt2pXZs2dz6tQpkz+l\nEI7l0KFD1K5dmy1btvDdd9+xaNGiTF0EQApBmqxcuZLKlStTo0YNWrVq9WCoyv3799O6dWvAOI9Q\nrlw5Dh48yKFDh2jevDkHDx6katWqhIWFUahQIZRSHDp0iKZNmwJw/PhxqlWrRnR0NLlz5zbzIwrh\nUJYsWfKgq4idO3cyYMCADNVVRGplikNDZnnSUJUnTpxgyJAhhIaGUrhwYQCyZcvG2LFjAWjevDnD\nhg3DxcWFihUrAnD27FkqVaoEQKdOnfDy8iJ37tx89NFHNvgkQji2mJgYBg4cyLfffkvDhg1Zvnw5\nTz/9tNmxbEYKgRXcLxCFChVizJgxDBo0iJ49ez74w/L09KRVq1YPnbtYuHDhg/uvv/46r7/+um1D\nC+GgLl26xLBhwzhx4gQjRoxg0qRJGbariNRyrE9rAhcXF2bNmmV2DCHEYwQEBNC5c2fCwsJYsWIF\nHTp0MDuSKeQcgRDC4Wit+fLLL2nSpAkFChTgu+++c9giAFIIhBAOJiwsjM6dO/P+++/zxhtvsG/f\nPkqVKmV2LFNJIRBCOIxTp07xwgsvsGrVKj7//HNWrFhBvnz5zI5lOjlHIIRwCKtXr6ZXr17kzJmT\nLVu20KiRDI1yn+wRCCEytdjYWEaOHEn79u2pXLkyQUFBUgQeYc0xi3MqpfYppQ4rpY4rpT6xTC+j\nlNqrlDqrlFqulMpurQxCCMd2/fp1mjdvztSpU+nfvz87d+6kZMmSZseyO9bcI4gCGmmtawA1gRZK\nqReBqcDXWusKwE2grxUzCCEc1L59+6hVqxZ79uxh4cKFfP/99+TIkcPsWHbJaoVAG8ItD50sNw00\nAlZapi8G3rBWBiGE49FaM3fuXF5++WWyZcvGb7/9Ru/evc2OZdeseo5AKZVVKXUIuAZsAc4Bt7TW\nsZaXhADFrZnB2ubMmUORIkWoWbMmZcuW5YcffjA7khAOKyIigr59+9K/f3/c3d0JDAykVq1aZsey\ne8oWfd8rpfIDa4CxwCKtdXnL9JLARq31f4b1Ukp5AV4Arq6ubr6+vg897+LiQvny5ZNcdlxcHFlT\nMzJNMg0fPpyqVavSt29fDh06RNu2bfnrr7+snis4OJjbt2+nev4nCQ8Px9nZOd3fN60kV8o4Yq6r\nV68yduxYzp49S/fu3enZs2ey/4/Z6/qCtGVzd3cP0lrXTvKFWmub3IBxwAggFMhmmVYP2JzUvG5u\nbvpRJ06c+M+0x7lz506yXpdaL7/8sg4ICNBaa/33339rFxeXZM2X1lzJ/fwp5e/vb5X3TSvJlTKO\nlsvPz08/9dRT2sXFRa9bty7F89vr+tKxsWnKBgTqZGyfrdlqqLBlTwClVC6gCXAS8AfuX8vdE1hr\nrQy2cPToUZ599lm01sycOZNWrVqZHUkIhxEfH89nn31Gy5YtKV68OIGBgQ+6gM/QYmJg/nwoX57c\nf/xh9cVZ84KyosBipVRWjHMRP2mtNyilTgC+SqnPgIPAgjQvKZFBi3PFxWGtQYsvXrxIeHg4zZs3\nx8nJibp16zJ79uyUL0sIkWK3bt2iR48erF+/nq5duzJ37lzy5Mljdqy0iY8HX18YNw6Cg6FuXbLE\nxFh9sVYrBFrrI8Dzj5l+HqhrreXa0pEjR2jcuDF+fn4PTb927RotWrSgefPmnD59mpUrV5Ili1y7\nJ0R6OXr0KB4eHvz555/MnDmTQYMGZewBZLSGdetg9Gg4dgyqV4e1a6F1a8J37LD64jNHFxOJ/HKP\nCAuz2pjFR48epUaNGv+Zvn//fjw9PRkxYgQDBgzgxo0bDwaoEUKkjY+PD/369cPFxYWAgAAaNGhg\ndqTU0xq2bjUKwL59UKEC+PhAp05gwx+P8jM1DY4ePUr16tX/M33//v0PCsTt27elCAiRDmJiYhg6\ndChdu3bFzc2NAwcOZOwisGcPuLtDs2Zw9SosWAAnTkCXLjYtApBZ9ghMsmzZssdOP336NDdu3GDN\nmjX069fPxqmEyHyuXLlCp06d2L17N8OGDePzzz/HycnJ7Fipc+CAsQewaRO4usLMmeDlBSZe9SyF\nwAqyZcvGN998Y3YMITKF3bt307FjR+7cuYO3tzeenp5mR0qdEydg7FhYtQoKFIApU2DQILCDE9xy\naMgKli5danYEITK8+02y3d3dyZs3L3v37s2YReD8eejZE6pVg82bjWLwxx/w4Yd2UQRA9giEEHbo\n7t27eHl54e3tTZs2bViyZAkuLi5mx0qZmBiYMAEmT4Zs2eDdd42Nvx2eM5RCIISwK8HBwXh4eHDs\n2DEmTpzIyJEjM17z63Pn4M03Ye9e6NbNOAxU3H67VZNCIISwG+vXr6d79+5kzZoVPz8/mjVrZnak\nlNEaFi+GwYONvQBfX+jc2exUScpgZVYIkRnFxcUxZswY2rRpQ7ly5QgKCsp4ReDmTWOj37s3uLnB\n4cMZogiA7BEIIUz2zz//8Oabb+Ln50fv3r2ZPXs2uXLlMjtWygQEQPfuxvUAkyfDiBGp69rGJFII\nhBCmOXjwIB4eHly+fJk5c+bw1ltvZayuIqKjjVZAn38O5cvD779D7aR7fbY3cmhICGGKxYsXU79+\nfWJjY9m1axdeXl4ZqwicPg316sHUqdCvHxw8mCGLAEghEELYWFRUFG+//Ta9evWifv36HDhwgLp1\nM1A/lFrDvHlQqxb8+SesXg1z59rNNQGpIYeGhBA2c/36dV599VX27t3LBx98wMSJE8mWLQNthkJD\n4a234OefoUkTo4VQsWJmp0oz2SNIo8OHD/PKK69QpUoVsmTJglKKcePGmR1LCLsTEBCAl5cXx48f\nZ+XKlUydOjVjFYFffzW6h964Eb780rhKOBMUAcgkewTD/IZx6OrjB6ZJ7djANYvUZHqLxAemiYyM\npHPnzixZsoS6desyZswYIiMjGT9+fIqXJ0RmpbXmyy+/ZOTIkRQvXhw/Pz8qV65sdqzki4yEjz+G\nr7+GypWNQlCzptmp0lWmKARm2bp1K7Vq1XpwfLN69er4+fk99oTX5MmTuXHjBv369WPNmjVcuXKF\nd955h2effdbWsYWwmbCwMPr27cuKFSto3749ffr0yVhFID4e2rUDPz8YONBoHZQ7t9mp0l2mKASJ\n/XIPs+LANMeOHaNatWoPHh84cIBatWoRExPDuHHjuHfvHtHR0fTs2RMfHx/69+/P7du38fHxoWfP\nnlIERKZ26tQpPDw8OH36NNOmTeO9995jhw1G20pXkyYZRWD2bHjnHbPTWE2yzxEopQoopaoqpcoq\npZKcTylVUinlr5Q6qZQ6rpQaapk+Xil1SSl1yHJ7LS0fwEwFCxbkyJEjAJw5c4bVq1fTpUsX5s6d\nS0REBPnz5yc8PJyKFSvSsGFDBg4c+OC+l5eXyemFsJ7Vq1dTt25dQkND2bp1K++//37GahoK4O9v\njB3ctSu8/bbZaawq0T0CpZQLMBDwBLID14GcgKtS6n/At1pr/yfMHgu8p7U+oJTKCwQppbZYnvta\na/1FunwCE3l6erJu3Tqee+45ChUqhI+PDwULFuTgwYPMnj2bHJaBJvz9/R+MWHbo0KHHDm8pRGYQ\nGxvL6NGjmTp1KnXr1mXlypWULFnS7Fgpd/UqeHpCxYowZw5ktCKWQkkdGloJLAFe1lrfSviEUsoN\n6K6UKqu1XvDojFrrK8AVy/0wpdRJwH6730sFZ2dn1q9f/5/pbdu2pVevXpQsWZJGjRpx6tQpXn75\nZcBoZXT/vhCZyfXr1+nSpQvbt2+nf//+zJgx48GPoQwlLs4oAnfuGOMJOzubncjqlNba+gtRqjSw\nE3gOGA70Au4AgRh7DTcfM48X4AXg6urq5uvr+9DzLi4ulC9fPsllp7bVkLWlNVdwcDC3b99Ox0SG\n8PBwnO3wD19ypYytc508eZLx48dz8+ZN3n33XVq2bGkXuZIrYa7SCxdSeulSTn34IVdbtDA5WdrW\nmbu7e5DWOunLnbXWyboBBYC6wCv3b8mczxkIAjwsj12BrBjnJyYCC5N6Dzc3N/2oEydO/Gfa49y5\ncydZr7O1tOZK7udPKX9/f6u8b1pJrpSxZa65c+fq7Nmz69KlS+ugoKBEX2v368vPT2ultO7d29Q8\nCaVlnQGBOhnb6WS1GlJK9QOGAiWAQ8CLwO9AoyTmcwJWAcu01qsthefvBM/PAzYkJ4MQwr5ERkYy\naNAgFixYQPPmzVm2bBkFCxY0O1bqhYQYg8hUrQqzZpmdxqaS22poKFAH+Etr7Q48j3Hi+ImU0URg\nAXBSa/1VgulFE7ysHXAsRYmFEKb766+/eOmll1iwYAFjxozhl19+ydBFQMXGGucFIiJgxYpMea1A\nYpJ7HUGk1jpSKYVSKofW+pRSqlIS8zQAugNHlVL3L/v9GPBUStUENPAn0D81wYUQ5tiyZQuenp7E\nxsaybt06WrdubXakNCuzYAHs3g3LloEDXt+T3EIQopTKD/wMbFFK3QQuJzaD1no38Lg2VxtTFjHR\nZWS8tsnpQNvgBL8Qj4qPj2fq1KmMGjWKqlWrsnr1aipUqGB2rLTbsIFnfH1hwADjmgEHlKxCoLVu\nZ7k7XinlD7gAflZLlQw5c+bkxo0bFCxY0KGKgdaaGzdukDNnTrOjCAdy+/Ztevbsydq1a/H09GTe\nvHnkycDdLj/w11/QowdhFSqQ9+uvzU5jmqQuKMsJDADKA0eBBVpru7hGvESJEoSEhHD9eqKnKoiM\njLTLjWZacuXMmZMSJUqkcyIhHu/YsWN4eHjwxx9/MGPGDAYPHpw5fnxFR0OnThAXx4lx43jBDrcT\ntpLUHsFiIAbYBbQEqmCcODadk5MTZcqUSfJ1AQEBPP/88zZIlDL2mkuIhHx9fenbty/58uXD39+f\nl156yexI6eeDD2DfPli5kogMfKI7PSTVaqiK1rqb1noO0AGQS2KFcAAxMTG8++67eHp6UqtWLQ4c\nOJC5isDq1TBjBgwZAu3bm53GdEkVgpj7d7TWsVbOIoSwA1evXqVx48ZMnz6doUOHsn37dooWLZr0\njBnFuXPQuzfUrQvTppmdxi4kdWiohlLqjuW+AnJZHitAa63zWTWdEMKm9uzZQ8eOHbl9+zbLli2j\na2ZrRRMZaZwXyJIFli+H7NnNTmQXEi0EWmv766RHCJHutNbMmjWL4cOHU7p0aTZv3vzQWBuZxvDh\ncOAArF0LpUubncZuJNVq6KnEntda/5O+cYQQtnb37l369+/PsmXLaN26NUuWLCF//vxmx0p/Pj7w\n3XcwYgS0aWN2GruS1KGhUCAEY2wBePgCMQ2UtUYoIYRtBAcH4+HhwbFjx/jss8/46KOPyJIl2eNV\nZRynT4OXF9SvDxMnmp3G7iRVCL4BGgJ7AB9gt5bLWoXIFDZs2EC3bt3ImjUrmzZtonnz5mZHso6b\nN8HDA3LkMM4LODmZncjuJFr6tdZDgZrACox+gw4qpT5XSiXdgF8IYZfi4uIYO3YsrVu3pmzZsgQF\nBWXeIhARYRwGOnsWfvoJ5ELMx0qyiwnLHoC/Uuog0AWYAJwF5lk5mxAind24cYNu3brh5+dH7969\nmT17Nrly5TI7lnXExkKXLrBnD/j6QqNEe813aEmdLM4DtAU6A4WB1UAtrfVFG2QTQqSjwMBAOnTo\nwJUrV/j+++/x8vLKHF1FPI7WxoDz69bBN98YTUbFEyW1R3AN49e/DxCMcYK4jlKqDsD9wWaEEPZt\n/vz5DBw4kCJFirB7927q1KljdiTrGjsW5s+H0aNh0CCz09i9pArBCoyN/7OWW0IaYw9BCGGnIiIi\nGDRoEAsXLqRZs2YsW7aMQoUKmR3Lur75Bj77DPr1g08/NTtNhpDUBWW9bJRDCJHOzp8/T4cOHTh4\n8CCjR49m/PjxZM2aya8RXb4chg6Ftm2NawYy66GvdJbUOYJugLfWOv4Jz5cDiloGoRFC2In//e9/\ntGtnDCOyfv16WrVqZXIiG9i6Fbp3hwYNjIvHsiV33C2R1JoqiNFkNAgIwhinOCfG+ASvYlxwNvJx\nMyqlSgJLgCJAPDBXaz3DcrXycqA0xlCVnbTWN9P8SYQQxMXF8cknnzBhwgRq1qzJqlWrKFvWAa77\nPHAA2rWDSpWME8SZtSWUlSR1aGiGUmoW0AhjDOLqQARwEuiutb6QyOyxwHta6wNKqbxAkFJqC9AL\n2Ka1nqKUGolRSD5M+0cRwrHduHGDrl278uuvv9KiRQtWr16deZuGJhQcDC1bwlNPgZ8fFChgdqJ0\nce0arFgBpUpZ/3Becq4jiAO2WG7JprW+Alyx3A9TSp0EimM0R21oedliIAApBEKkSWBgIO3bt+fq\n1avMmTOHChUqOEYRuHoVmjeHuDjYvBmKFzc7UZrcuQM//2wc2dqyxfhY48c/hbWP7NmkUxGlVGng\neWAv4GopEveLxdO2yCBEZqS1Zt68eTRo0ACA3bt3Z+7rAxK6fdvYE7h6FTZuhGcfbdiYMURGwpo1\n0LEjuLpCz55w8qTRN96RI/Dqq4kPx5selLW7DlJKOQM7gIla69VKqVta6/wJnr+ptf7PvpxSygvw\nAnB1dXXz9fVN1fLDw8NxdnZOXXgrklwpI7n+KyoqiunTp+Pn50ft2rUZPXo0Li4upudKTHrlyhId\nTbUPP8Tl6FGOTZrEP3Xr2kWu5IqLg0OHCrB169Ps2lWYu3ezkT9/NO7u12jU6BpVq9550OApLdnc\n3d2DtNa1k3yh1tpqN8AJ2AwMTzDtNEZLI4CiwOmk3sfNzU2nlr+/f6rntSbJlTKS62Hnzp3TNWvW\n1IAeM2aMjo2NtYtcSUmXXLGxWnfooDVovXRp2t9P22Z9xcdr/fvvWg8ZorWrqxE/b16te/bU2s9P\n65iY9M8GBOpkbKuT1b5KKeUKTAKKaa1bKqWqAPW01gsSmUcBC4CTWuuvEjy1DugJTLH8uzY5GYQQ\nhvXr19OjRw/A6EH09ddfNzmRDWltjDO8ciV8+SV062Z2oiQdP24c8/fxgfPnjU5QX38dunaF116z\njwZOyT1H8APGL/tilsdngGFJzNMAo8fSRkqpQ5bbaxgFoKlS6izQ1PJYCJGE2NhYRo4cSZs2bShT\npgxBQUGOVQTAuGL422+NA+jDh5ud5on++gumToUaNeC552DyZChXDhYtgr//hlWroH17+ygCkIxW\nQxaFtNY/KaU+AmMge6VUXGIzaOMisyedsWqcgoxCOLyrV6/SpUsXduzYgZeXFzNmzCBnzpxmx7Kt\nuXONPoR69IAp9vf78fp1o7mnt7fR4SnAiy/CjBlGn3dFipibLzHJLQR3lVIFMfoXQin1InDbaqmE\nEA/s2LGDLl26cPv2bRYvXvzgsJBDWbgQBgwwjqXMn28MPm8HHtfcs2pVYxC0Ll0go1zLl9xCMBzj\n2H45pdQejC6pO1gtlRACrTXTpk3j448/ply5cvz666+Zc0D5pMyYAcOGGdcLrFhh+ghjkZGwaZPx\ny3/DBuNxqVLG0aquXSEjfkVJFgKlVBaMbiVeBSphHO45rbWOsXI2IRzWrVu36NmzJ+vWraNjx47M\nnz+ffPnymR3LtrQ2flqPGWMMNentbZxpNUFcHPj7GxFWrzYuYShcGPr2NTb+9epl7P7tknNlcbxS\n6kutdT3guA0yCeHQDhw4QIcOHbh48SIzZsxg8ODBjnGBWEJaw4cfwrRpRkdyCxfavBM5rWHfPmPj\n/9NPxnVrefMaXRp17QqNG2eefu2S+zF+VUq1B1Zb2qYKIdKZ1pr58+czePBgChcuzM6dO6lXr57Z\nsWwvPh4GDoTvvzdGGZs1y6bnBE6cMDb+95t7Zs/+b3PP11+3n5Y+6Skl5wjyALFKqUiMw0Naa+1g\n+6pCWMe9e/d4++23WbJkieMMIPM4sbHQuzf8+KOxRzB5sk2OuVy9moOpU40CcOSIUXcaNzYGOGvX\nDvLnT/o9MrJkFQKtdV5rBxHCUZ05c4b27dtz/Phxxo8fz+jRozP/ADKPExUFnp5GxzsTJ8LHH1t1\ncQ839zT2vF58EWbONPr9sefmnuktuVcWv/K46VrrnekbRwjHsmLFCvr27Uv27Nnx8/OjWbNmZkcy\nx927xk/vLVuMVkJDhlhlMWFhRnNPb++Hm3v27Xuejz8um2Gae6a35B4aGpHgfk6gLsZANY3SPZEQ\nDiA6OpoPPviAGTNm8OKLL/KwV2kXAAAgAElEQVTTTz9RsmRJs2OZ4/ZtaNUKfvsNFiyAPn3S9e2j\nov5t7rl+/eObewYEXHCMAXyeILmHhlonfGwZfexzqyQSIpO7cOECnTp1Yu/evQwZMoRp06aRPXt2\ns2OZIzQUWrSAw4eNs7OdOqXL295v7unjY3TncL+5Z79+xtGnjN7cM72ltvFTCPBcegYRwhFs3LiR\n7t27ExMTw4oVK+jQwYGvy7x8GZo2NZrmrF1rXDWcBvebe/r4GGPY32/u6eFhbPwzU3PP9JbccwTf\nYOleAqOjuprAYWuFEiKziY2NZezYsUyePJkaNWqwYsUKKlSoYHYs8/z5p7FlvnbNOG7TsGGq3+rE\nCWPj7+1tv7172rvk1sfABPdjAR+t9R4r5BEi07ly5Qqenp7s2LGDfv36MXPmTMcYRvJJTp2CJk3g\n3j3YuhVeeCHFb3HhAvj6Ghv/w4f/be45ZoxxztkyPo9IpuSeI1h8/75SqgDgoGe1hEgZf39/PD09\nCQsLc9wO4xLIExxsnAdQCgICoHr1ZM97/boxDIG3N+zebUxz1Oae6S25h4YCgDaW1x8Criuldmit\n7bdDcCFMFB8fz6RJkxg3bhwVK1Zk27ZtVK1a1exY5tEafH2p+e67UKCAsSdQsWKSs91v7unjA7/+\nmnF797R3yT005KK1vqOU6gcs0lqPU0odsWYwITKq0NBQunXrxubNm+natStz5syxy/GDbSYoCIYO\nhT17iKxYEadffzXabz7B/eaePj6wbt2/zT3ffx/efDNj9u5p75JbCLIppYoCnYBRVswjRIb222+/\n0blzZ65du8b333+Pl5eX43UYd9+1a8bVwQsXGm03FywgqHRpGj6mCMTFGUeKvL0fbu6ZWXr3tHfJ\n7cnpU4yhKoO11vuVUmWBs4nNoJRaqJS6ppQ6lmDaeKXUpUeGrhQiw9Na89VXX/Hqq6+SPXt2fv/9\nd/r37++YRSA62hhPuEIFWLzYGFLyzBnjQrEEncfdb+45bBiUKGGcP16xAtq2BT8/o3XprFlQv74U\nAWtL7sniFcCKBI/PA+2TmO0HYBaw5JHpX2utv0hBRiHs2q1bt+jTpw9r1qyhXbt2LFy4kPyZvZey\nJ9m0Cd59F06fNtpufvUVVKr00EtOnvy3d89z5xyjd097l9yTxZ8DnwERgB9QAximtf7xSfNorXcq\npUqnQ0Yh7NaZM2fo27cvFy5c4KuvvmLYsGGOuRdw5oxRADZuNE4C//LLQxeI3W/uOW+eG8HB/zb3\nHDXKMXr3tHfJPUfQTGv9gVKqHcZVxR0Bf+CJhSARg5RSPTCuTXhPa30zFe8hhKm01syZM4chQ4bg\n6urquGMH3LkDEyYYHcXlygVffAGDB0P27I9t7lm5ss4Qg7k7GpWccWaUUse11lWVUvOAVVprP6XU\nYa11jSTmKw1s0Fo/Z3nsCoRiXKU8ASiqtX5sD1NKKS/AC8DV1dXN19c3+Z8qgfDwcLtssSG5Usae\nct27d48vv/yS7du3U6tWLcaOHYuLnV3BZPX1FR9PET8/ys6fj9OtW1xt2ZLzfftyO2dh9uwpxLZt\nT7N//1PExytKlbpLkyZ/06jRNfLlC7Wb7zEhe/r7elRasrm7uwdprWsn+UKtdZI3YApwCjgIOGEM\nXr83GfOVBo6l9LlHb25ubjq1/P39Uz2vNUmulLGXXEePHtWVKlXSWbJk0RMnTtTbtm0zO9JjWXV9\n7dmjtZub1qB1/fo66rdAvWaN1p06aZ0rlzH5mWe0/vBDrQ8f1jo+3ka50sBec2mdtmxAoE7GNja5\nJ4tHKqWmAne01nFKqXtA2xQWJ5RSRbXWVywP2wHHEnu9EPZk8eLFvP322+TLl4+tW7fi7u5OQECA\n2bFs4+ZNY9R2Hx/Ytg1dvDjHP1rG11c9WdVSPWju2afPv7172nB0SZFGyT1ZnBsYCDyDcbimGFAJ\n2JDIPD5AQ6CQUioEGAc0VErVxDg09CfQPw3ZhbCJe/fuMXjwYBYuXEjDhg3x8fGhiCMc4L571+jA\n38fHaA0UE0NkyQpsq/cZg88N44/JeXB2Nnr3zGyDuTua5H5tizAGoqlveRyC0Zz0iYVAa+35mMkL\nUpROCJOdPn2ajh07cuzYMUaPHs348eMz9zCS0dGwebOx8V+7Fu7dI8a1OHvdhjAtxJN1F2uR/W/F\n66/D59LcM9NIbiEop7XurJTyBNBaRyiHbCMnHMny5cvp168fOXPmZNOmTTRv3tzsSNYRFwc7dvw7\nisvNm8QVKMjRat35JtSTRedeRl3PQqNGsOATYw9AmntmLsktBNFKqVxYxiRQSpUDoqyWSggTRUVF\nMXz4cL799lvq16/P8uXLKVGihNmx0pfWsH//v6O4XLlCfB5nzlZ9g/nhnkw/0ZTYvU688AJMHyLN\nPTO75BaCcRgXkpVUSi0DGgC9rBVKCLOcP3+eTp06ERQUxPvvv8+kSZNwcnIyO1b6+ecf+PZb+OEH\nOHcOnT07F6u9xjJXTyYfbUXYvtxUrgzjPzN69yxXzuzAwhaSLASWQ0CnAA/gRUABQ7XWoVbOJoRN\n/fzzz/Tq1QulFD///DNt26a4YZz9unTJ6O5hzhy4e5fr1Ruzss4oPj3ajqtB+XnmGXjnfaPFT/Xq\n0rePo0myEGittVLqZ621G/CLDTIJYVMxMTGMHDmSr776itq1a/PTTz9RpkwZs2OljzNn4PPP0UuW\noOPi2VemC+9d+5DfjlSjUCHoZGnuWb++NPd0ZMk9NPQ/pVQdrfV+q6YRwsYuXrxI586d+f333xk0\naBBffPEFOXLkMDtWmjmfPo2ePRtWrSI2S3aW5XiLT2LeJ/TvMrRrB6M9jd4+M9NRL5F6yS0E7sAA\npdSfwF2Mw0Naa538ceaEsDMbN26ke/fuxMTEsHz5cjp16mR2pLTRGgICCB89mdq/bSEsSz6+0SP5\nLstQajdzfdDcM3dus4MKe5PcQtDSqimEsKGYmBjGjBnD1KlTqVGjBj/99BMVkzFsot2Kj+f6gnXE\nTJhCsYt7CceVz9RkTr78Nm17uHBUmnuKJCRaCJRSOYEBQHngKLBAax1ri2BCWENISAhdunRhz549\n9O/fn6+//ppcGfSKqNArMRwZ6U3ZFVMpHXGS85Th89LfkvudXtQpe5Ap7e2rIzxhv5LaI1gMxAC7\nMPYKqgBDrR1KCGvYtGkT3bt3JyoqCm9vbzw9H3fxu30LC4P1KyK58+VcXjvxBY24yOkc1fi54zKq\nTejEB5WM/9IBAdEmJxUZSVKFoIrWuhqAUmoBsM/6kYRIX7GxsYwdO5bJkydTrVo1VqxYQaVHRs2y\nZ1FRxtCN9wdzj4+AC1kmcbdEBf748DsqvvMalbJIe0+RekkVgpj7d7TWsdKrhMhoLl26hKenJ7t2\n7eKtt95ixowZGeJQUMJeH1auhFu3oFAh6N0bunbNSaHSB3m6eFGzY4pMIqlCUEMpdcdyXwG5LI/v\ntxrKZ9V0QqTB5s2b6datGxEREfz444+8+eabZkdKlNYQGGiM6GXp9eFB756enkbvnv8295QiINJP\nooVAa52Ju1kUmVVsbCzjx49n0qRJVK1alRUrVvDss8+aHeuJTp40fvn7+EBw8L+DuXt6QqtW0run\nsD7pPVxkKleuXMHT05MdO3bQt29fZs6cSW47bDh/8aIxmLuPDxw8aFzV26gRfPSR9O4pbE8Kgcg0\ntmzZQrdu3QgPD2fx4sX06NHD7EgPCQ39dzD3XbuMaS+8ANOnG717FpWjPcIkUghEhhcXF8cnn3zC\nZ599RuXKlfH396dKlSpmxwIgPNwY38XbG379FWJjoXJlmDDBOPQjvXsKeyCFQGRoV65c4c0338Tf\n359evXoxa9Ys8uTJY2qmqChjkC9vb6O5Z0QEPPMMvPee9O4p7JPVCoFSaiHQCrimtX7OMu0pYDlQ\nGmPM4k5a65vWyiAyt4SHghYtWkSvXr1MyxIXBzt3Ghv/R5t7Su+ewt5Z80/zB6DFI9NGAtu01hWA\nbZbHQqRIbGwso0ePpnnz5hQuXJj9+/ebUgS0htOn8zJ8OJQsaZzs9fWF1q1h40a4fBlmz4aXXpIi\nIOyb1fYItNY7lVKlH5ncFmhoub8YCAA+tFYGkfmEhITQtWtXdu3aZVqroFOnjNY+3t4QHOxG9uzw\n2mvQVXr3/I+wqDAOXj1I4OVAAi8HcuCvA1S/Xp3KhSpTpXAVqhSuQsWCFcmRLeN3/Z2RKa219d7c\nKAQbEhwauqW1zp/g+Zta6wJPmNcL8AJwdXV18/X1TVWG8PBwnJ2dUzWvNUmulAkPD+fYsWNMnjyZ\n6Ohohg8fTtOmTW22/GvXcrB9+9Ns3/40Z8/mJUsWTc2at2jQ4ALNmoXh7GxffTGa8T1GxkUSHB7M\n6bDTnA4/zZmwM1y4dwFtDHWOaw5XimQvQmhsKJcjLj+YnoUsFM9VnGdyP0Op3KUonac0pXKX4pnc\nz5Aza06bZLfXv3tIWzZ3d/cgrXXtpF5nt4Ugodq1a+vAwMBUZQgICKBhw4apmteaJFfyxcTE0L17\nd5YvX27TbqMf19yzbl3jl//95p72uL7A+rmiYqM48veRB7/0A68EcvzaceJ0HABFnItQp1gdaher\nTZ1idXAr5sbTeZ5+kCsiJoIzN85wMvQkJ66feHA7+89ZYuONoqpQlMpfythzKGTsPTQu25hnXJ5J\n989jr98jpC2bUipZhcDWrYb+VkoV1VpfUUoVBa7ZePkig7lw4QJdunTh999/Z8CAAXz11VdW7Sso\nseaeXbpA+fJWW7Tdu373OnOC5rDm1BqO/n2UmHijK7JCuQtRp1gd2lZqS+1italdrDbF8hZL9L1y\nOeWiRpEa1ChS46HpMXExBP8TzInrJx4qEtvObyMqLgoAt6JueFT2wKOyB88Wst8rxjMSWxeCdUBP\nYIrl37U2Xr7IQNatW0evXr0e9B76ySefWGU5j/buGRFhnPwdPtz49e/ozT2P/H2EGf+bwbKjy4iK\ni+LlZ17mvXrvPdjoP+PyDOnVIaVTVicqF65M5cKVH5oeFx/HmRtnWH9mPatPrmbU9lGM2j6KyoUq\nPygKzxd5Pt1yOBprNh/1wTgxXEgpFQKMwygAPyml+gIXgI7WWr7IuKKjo/nwww+ZPn06tWrVYvny\n5YSEhKTrMp7Uu2evXsbG39Gbe8bFx7HhzAZm7J2B/5/+5HbKTe+avRnywpD/bKRtIWuWrA8KxAcN\nPuDSnUv8fOpnVp9azZTdU5i4ayKlXErR7tl2eFT2oH7J+mTNIl2lJZc1Ww09adSPxtZapsj4/vjj\nDzp37sz+/fsZPHgw06ZNI0eOHOlSCJ7Uu+cbbxgbfxnMHe5E3WHhwYV8s+8bzt88T8l8JZnaZCr9\navXjqVxPmR3vgeL5ijOw7kAG1h1I6L1Q1p9ez+pTq/ku8Dum753O03me5o1Kb9CucjsalWlE9qzZ\nzY5s1+TKYmE3Vq9eTZ8+fQBYtWoVHh4e6fK+j+vd87XX/u3dU5p7QvA/wXyz9xsWHVpEWHQYDUo2\nYErjKbSr3I5sWex7M1EodyF6P9+b3s/3JiwqjE3Bm1h9cjXex7yZe2AuLjlcaFWxFf1q9aNh6YZm\nx7VL9v0NC4cQERHBBx98wKxZs6hTpw7Lly+nTJkyaXrP+717envDoUPGMf5GjWDkSKN3zwJJtlXL\n/LTWbP9jOzP2zmDDmQ1ky5KNzs91ZugLQ6ldLMmGJnYpb468dKraiU5VOxEZG8m289tYfXI1a0+v\nZdnRZbQo34Ipjaf85yS1o5NCIEx14MABunXrxsmTJ3n33XeZMmUK2bOnbjf+Sc09v/7aaO5ZLPGG\nLA4jIiYC76PeTN87nWPXjlE4d2FGvzKat2u/TdG8macL1JzZcvJ6xdd5veLrzIqZxez9s5m0axLP\nz3meN6u/yQT3CZTOX9rsmHZBCoEwRWxsLFOnTmX8+PE8/fTTbN68mWbNmqX4fR7X3PPZZ+HTT41D\nP47c3DOh0Huh/HLmF9adWcfm4M3cjblLDdcaLGyzEM9qnuTMZpsLt8ySyykX79d/n77P92XqnqnM\n2DuDn47/xDu132HUK6PMjvcfkbGR7L6wm63nt+IW52b15UkhEDYXHBxMjx49+P333+ncuTPffvst\nTz2V/BOR0dFGc8+EvXveb+7p6Qk1ajh2c8/7zt44y9rTa1l3eh17Lu4hXsdTPG9xetToQeeqnXml\n1CsO19yyQK4CTGkyhUF1BzE+YDwz981kwcEFdCrWiTrRdciT3Zyea+N1PIevHmbL+S1sOb+F3Rd2\nExkbiVMWJz6pYp1m0wlJIRA2o7Vm3rx5DB8+HCcnJ7y9vfH0fFLjsofFxYG/v7HxX7UKbt6EggWN\n5p6entCggWM39wSjyefeS3tZd3odPgd9uLDjAgDVXasz6uVRtK3UllpFazncxv9xSuQrwfw283mv\n3nt8vP1jFpxawC/f/ML4V8fT5/k+OGW1fvOxv279xZbzW9h6fivb/thG6L1QAJ57+jkGuA2gabmm\nvFLqFQJ/S12vCikhhUDYxNWrV+nXrx+//PILjRs3ZtGiRZQsWTLRebSGoCBj4790aT1CQyFPHmjX\nTpp73ncv5h5bz29l7am1bDi7gWt3r5EtSzaq56vOe6+8R5tKbeQ4eCIqF67Mms5rmLV2Fr43fBnw\nywC++t9XTGw0kfaV26dr0bwVeQv/P/wfbPzP/nMWgKLORXmtwms0KdOEJmWbmHKeRgqBsLo1a9bg\n5eVFeHg406dPZ/DgwWRJ5Of7w717Gs0969QJY/bsHNLcE6Ofnx+P/Mja02vZcn4LkbGR5MuRj9cq\nvEabim1oWaElh/53iIYvNDQ7aobxnMtz7Gqziw1nNjBy20g6ruhI3eJ1mdpkarKbnMbFx3Er8hah\n90K5EXGDG/duEHovlHM3z7H1/Fb2X95PvI7HObszDUs3ZGCdgTQt15TKhSqbvpcmhUBYzZ07dxg6\ndCg//PADtWrVYunSpU8cQjIk5N/mngcP/re55+HDx+y2UzBby5olKyO2jCBvjrz0e74fbZ9tyyul\nXpGLptJIKUXrSq15rcJrLDm8hLEBY3Ff7E7L8i3pUaMHd6LucOPeDW5E3PjPxv5GxA1uRtx80KNq\nQllVVuoWr8vol0fTpGwTXijxgt19V1IIhFXs3LmTHj16cPHiRUaNGsXYsWP/0yz0xo2Hm3tqbTT3\nlMHcE5ctSzYODzhMiXwlTP8lmRllzZKV3s/3pstzXR40Od0UvOnB87mdclMwV0EK5S5EwdwFecbl\nmYceP3rf1dmV3E72vRsrhUCkq6ioKMaMGcMXX3xB2bJl2bVrF/Xr13/wfHi40dLH29sY1/d+756f\nfCLNPVOipEvi51dE2t1vcurl5sUfN/94sGHP5WS93m/NIoVApJsjR47QvXt3jhw5gpeXF19++SXO\nzs5ERz88mPu9e9K7p8g48uXIl+mvRJZCINLs0qVLTJ06lTlz5pA/f37Wr19Py5at2LXr38Hc7zf3\n7NlTevcUwt5IIRCpFhISwpQpU5g3bx7x8fF0794DT88v8PMrQP/+xuDt0txTCPsnhUCk2MWLF5ky\nZQrz588nPj6eN974kGLF3mPTpgIsWmQ092zZ0tj4S3NPIeyfFAKRbBcvXmTy5MksWLCAuLii1K79\nI3fvtmHlyhwoBe7u8OGH0runEBmNFAKRpAsXLjB58mTmz19DfLwHrq5HuXq1Anv3KurUkd49hcjo\nTCkESqk/gTAgDojVWmfMzs8zub/++otPP/2KxYtvEh/fGaVmER+fFRcXePttae4pRGZh5h6Bu9Y6\n1MTliyc4c+ZPBg/+hS1bCqH1JCAPxYrF8uabWenaVXr3FCKzkUNDAjB69/zhh/N88kkkFy/mAwaS\nM2c4HTpovLygQYNs0txTiEzKrEKggV+VUhqYo7Wea1IOh6Y1rFt3iWnTQti7tzSxsWWBp6lY8SSj\nRsXj6VlImnsK4QCU1v/tJMnqC1WqmNb6slLqaWALMFhrvfOR13gBXgCurq5uvr6+qVpWeHg4zs7O\naY2c7szMdeRIND/+GMfhw1WIji4FRJM37x4aNPiLDh1yUK6c/XXyI99jykiulLHXXJC2bO7u7kHJ\nOgertTb1BowH3k/sNW5ubjq1/P39Uz2vNdk61/Hjt3WnTvt03rxntLEvEKednf+nPTw26sOHL5iW\nK7kkV8pIrpSx11xapy0bEKiTsR22+aEhpVQeIIvWOsxyvxnwqa1zOIKQkAgmTDjB6tU5CQ2tDNQh\ne/bDNGmykY8/Loe7+wtmRxRC2AEzzhG4Amss3edmA7y11n4m5MiUbt2KZcqUE/z4YzyXLlUF3Mia\n9SwvvriJ4cOL06FDDZTK3B1oCSFSxuaFQGt9HpAtUTqKitLMmHGK+fPvcvZsZaA6SoVQrdpW3nnn\nKfr2dcPJqYLZMYUQdkqaj2ZQsbGaxYvP8803Nzh6tCLx8ZWBUMqW3UPv3jkZPrwuuXO3NDumECID\nkEKQgRjNPUMszT3LEBtbDnClaNF9dO2q+eijOhQs2MzsmEKIDEYKQQawe/d1Jk48T0BAUSIjnwGe\n5qmn9vLGG4cZM+Z5SpduZHZEIUQGJoXATh0/fptPPz3Dpk0uhIVVBAri7Lyf9u1PMHZsVapXf9ns\niEKITEIKgR0JCYng009PsGZNDkJDqwB1yJHjME2bbuTjj8vTsKE09xRCpD8pBCa7eTOGyZOP4+2t\nuXTpOf5t7rmZ998vhodHdWnuKYSwKikEJoiMjMfbO4y33tpPcHAVoCZKXaJatW0MGvQUffq4kS2b\nNPcUQtiGFAIbiYvTLFwYzKxZ/3DsWCXi41sDoZQv/xt9+uRi2LA65MrVwuyYQggHJIXAirSGn3++\nwLRpl9i/vwyxsRWAcIoX389LL13ku+88KFCgqdkxhRAOTnqYt4KdO/+mefPfyZXrLzw8nuH3391w\ncQmmX7+tXLgQRUiIOwMGPEOBAvbZ26EQwrHIHkE6OXr0Jp9+ego/v6cID68EFCZv3iDatDnNuHHP\nUbXqS2ZHFEKIx5JCkAZ//RXOp58e5+efc/PPP1WBeuTMeYTmzTczalR5Xn65jtkRhRAiSVIIUuif\nf6KYOPEovr6Ky5erAy+QLVswDRpsYcSIkrRpUw2lqpsdUwghkk0KQTJERMTx5ZdHWbQokvPnqwG1\nyZLlEjVrBjB4cCF69qxB1qzlzY4phBCpIoXgCeLiNHPmnOC7725x/HhltK4J/EPFinvp1y8PQ4bU\nIkeO4mbHFEKINJNCkIDW8NNPZ/nqq6sEBZUnLq4qEE7Jkgfo1i0rH35YCxcX6eBNCJG5SCEAtm27\nyKRJf7B7d0mioysApShcOIj27c8wZkxNihV7xeyIQghhNaYUAqVUC2AGkBWYr7WeYusMBw78zaef\nnmbr1sLcvVsZKI6Ly0E6dPiTsWOfo1KleraOJIQQpjBj8PqswGygKRAC7FdKrdNan7D2ss+fv8X4\n8UdZvz4vt25VB1zJles4rVptZ8yYStSt62btCEIIYXfM2COoCwRbxi5GKeULtAWsUghu3oxm8ODf\nWLEiG3//XRN4GSen87z6agAjR5aiRYuqQFVrLFoIITIEMwpBceBigschgFU62m/cOIDt2xsAecia\n9TJubnt4911XPD0rkyVLWWssUgghMhyltbbtApXqCDTXWvezPO4O1NVaD37kdV6AF4Crq6ubr69v\nipc1bdp1jh7NTdu292jbtiDZstlP10rh4eE4O9tfX0OSK2UkV8pIrpRLSzZ3d/cgrXXtJF+otbbp\nDagHbE7w+CPgo8TmcXNz06nl7++f6nmtSXKljORKGcmVMvaaS+u0ZQMCdTK2y2b8RN4PVFBKlVFK\nZQe6AOtMyCGEEAITzhForWOVUoOAzRjNRxdqrY/bOocQQgiDKdcRaK03AhvNWLYQQoiH2c/ZUyGE\nEKaQQiCEEA5OCoEQQjg4KQRCCOHgpBAIIYSDs/mVxamhlLoO/JXK2QsBoekYJ71IrpSRXCkjuVLG\nXnNB2rKV0loXTupFGaIQpIVSKlAn5xJrG5NcKSO5UkZypYy95gLbZJNDQ0II4eCkEAghhINzhEIw\n1+wATyC5UkZypYzkShl7zQU2yJbpzxEIIYRInCPsEQghhEhEpi4ESqkWSqnTSqlgpdRIGy63pFLK\nXyl1Uil1XCk11DJ9vFLqklLqkOX2WoJ5PrLkPK2Uam7lfH8qpY5aMgRapj2llNqilDpr+beAZbpS\nSs20ZDuilKplpUyVEqyXQ0qpO0qpYWasM6XUQqXUNaXUsQTTUrx+lFI9La8/q5TqaaVc05RSpyzL\nXqOUym+ZXlopFZFgvX2fYB43y/cfbMmurJArxd9bev9/fUKu5Qky/amUOmSZbsv19aTtg3l/Y8kZ\ntCAj3jC6uD4HlAWyA4eBKjZadlGgluV+XuAMUAUYD7z/mNdXseTLAZSx5M5qxXx/AoUemfY5MNJy\nfyQw1XL/NWAToIAXgb02+u6uAqXMWGfAK0At4Fhq1w/wFHDe8m8By/0CVsjVDMhmuT81Qa7SCV/3\nyPvswxggSlmyt7RCrhR9b9b4//q4XI88/yUw1oT19aTtg2l/Y5l5j6AuEKy1Pq+1jgZ8gba2WLDW\n+orW+oDlfhhwEmOs5idpC/hqraO01n8AwRj5baktsNhyfzHwRoLpS7Thf0B+pVRRK2dpDJzTWid2\nEaHV1pnWeifwz2OWl5L10xzYorX+R2t9E9gCtEjvXFrrX7XWsZaH/wNKJPYelmz5tNa/a2NrsiTB\nZ0m3XIl40veW7v9fE8tl+VXfCfBJ7D2stL6etH0w7W8sMxeC4sDFBI9DSHxjbBVKqdLA88Bey6RB\nlt27hfd3/bB9Vg38qpQKUsbY0ACuWusrYPyhAk+blA2MUesS/ge1h3WW0vVjxnrrg/HL8b4ySqmD\nSqkdSqmXLdOKW7LYIvLCa8wAAAVESURBVFdKvjdbr6+Xgb+11mcTTLP5+npk+2Da31hmLgSPO45n\n0yZSSilnYBUwTGt9B/gOKAfUBK5g7JqC7bM20FrXAloCA5VSryTyWptmU8bwpW2AFZZJ9rLOnuRJ\nOWy93kYBscAyy6QrwDNa6+eB4YC3UiqfDXOl9Huz9ffpycM/Nmy+vh6zfXjiS5+QId2yZeZCEAKU\nTPC4BHDZVgtXSjlhfMnLtNarAbTWf2ut47TW8cA8/j2UYdOsWuvLln+vAWssOf6+f8jH8u81M7Jh\nFKcDWuu/LRntYp2R8vVjs3yWk4StgDcthy+wHHq5YbkfhHH8vaIlV8LDR1bJlYrvzZbrKxvgASxP\nkNem6+tx2wdM/BvLzIVgP1BBKVXG8iuzC7DOFgu2HH9cAJzUWn+VYHrCY+vtgPutGdYBXZRSOZRS\nZYAKGCeorJEtj1Iq7/37GCcbj1ky3G910BNYmyBbD0vLhReB2/d3X63koV9q9rDOEiwvJetnM9BM\nKVXAclikmWVaulJKtQA+BNpore8lmF5YKZXVcr8sxvo5b8kWppR60fJ32iPBZ0nPXCn93mz5/7UJ\ncEpr/eCQjy3X15O2D5j5N5aWs9/2fsM4234Go7qPsuFyX8LYRTsCHLLcXgOWAkct09cBRRPMM8qS\n8zRpbJWQRLayGC0yDgPH768XoCCwDThr+fcpy3QFzLZkOwrUtmK23MANwCXBNJuvM4xCdAWIwfjV\n1Tc16wfjmH2w5dbbSrmCMY4T3/87+97y2vaW7/cwcABoneB9amNsmM8Bs7BcWJrOuVL8vaX3/9fH\n5bJM/wEY8Mhrbbm+nrR9MO1vTK4sFkIIB5eZDw0JIYRIBikEQgjh4KQQCCGEg5NCIIQQDk4KgRBC\nODgpBCJTU0rFKaM3yeNKqcNKqeFKqVT/3SulPk5wv7RK0LNlEvMNU0r1SMFysiuldloufhLCqqQQ\niMwuQmtdU2tdFWiK0V57XBre7+OkX/Iwy8a8D+Cd3Hm00fHaNqDz/9u7f9cooiiK49+TIJp6UbCw\nUyst10aRaGFlZWtjoyJipUjsLdLZKAS1s09SKbESEYRIQvBXsNE/II1WmsTkWtw3JBmzmAm7KZzz\nqYZh3rxh2Z277w1zXtP+zJpyIbDWiIzUuEaGoUnSsDLP/10JR7sOIGm0/BufkvRZ0oSkIUnjwEgZ\nYVSZPsOSnpQRx0tJI9t0fZ6Mzfhdzv9K0oPSx6KkrqRJZab8/U3tpoHLg/tEzJILgbVKRHwlv/eH\nyDdgf0REF+gCV0vsAWQ2zm3gJBmedikixtgYYVQ36GPAozLi+E6+oVp3Gpir7VuJiLPABBklcBM4\nAVyR1CnHfCzXZTZQLgTWRlVq4wUyw2WBjAHukDd2gNnIbPw1MqrgTI9zfYuIhbI9Ry5wUncYWKrt\nq3J0PgCfIjPql8nFRY4AlL5Xqmwos0HxgyhrlRIotkYmOwq4FREztWNG+TvOt1cWy/Km7TVgu6mh\nn8CBHu3Wa+dYZ+vvcj/wq0ffZn3hEYG1hqSD5FTMw8iQrRngRokERtLxksgKcKokYQ6RD2zflP2r\n1fENLAJHd3G9HWApIlabtjVrwiMC+9+NlKmffeTCLc+AKvr3KTmVM1+igZfYWB7wLTBOPiN4Ta7b\nAPAYeC9pnkzR3IkXpd+mzgHPd9HOrBGnj5rVlKmhOxFxsY/nnALuxtalEf/VZhK4FxFf+nUdZtvx\n1JDZ3hgjHxrvSFmcZdpFwPaCRwRmZi3nEYGZWcu5EJiZtZwLgZlZy7kQmJm1nAuBmVnLuRCYmbXc\nHwI2WCKzH16KAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10d935358>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy import integrate\n",
    "import scipy.optimize\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "#Input data\n",
    "depth = np.linspace(0, 2000, num=21)\n",
    "density = np.array([1025, 1026, 1026, 1030, 1030, 1031, 1900, 2190, \n",
    "                    2200, 2230, 2235, 2240, 2275, 2305, 2310, 2308, \n",
    "                    2310, 2305, 2310, 2323, 2319])\n",
    "density_brine = 1031. * np.ones_like(density)\n",
    "porosity = np.array([0.305, 0.297, 0.286, 0.281, 0.285, 0.293, 0.307, 0.305, 0.298])\n",
    "\n",
    "#Compute S_v\n",
    "Sv = integrate.cumtrapz(density * 9.81, depth, initial=0)\n",
    "#Compute S_hydro\n",
    "hydro = integrate.cumtrapz(density_brine * 9.81, depth, initial=0)\n",
    "#Compute P_p\n",
    "phi_0 = 0.38\n",
    "beta = 3e-2 / (1e6)\n",
    "Pp = Sv[-len(porosity):] + (1. / beta) * np.log(porosity / phi_0)\n",
    "#Compute sigma_eff\n",
    "sigma_eff = Sv[-len(porosity):] - Pp\n",
    "#Create plot\n",
    "plt.plot(depth, Sv / 1e6, 'k', label=\"$S_v$\") \n",
    "plt.plot(depth, hydro / 1e6, 'b', label=\"$S_{hydro}$\") \n",
    "plt.plot(depth[-(len(porosity)):], Pp / 1e6, 'r', label=\"$P_p$\")\n",
    "plt.plot(depth[-(len(porosity)):], sigma_eff / 1e6, 'g', label=\"$\\sigma_{eff}$\");\n",
    "plt.grid()\n",
    "plt.legend(loc=2)\n",
    "plt.xlabel(\"Depth (m)\")\n",
    "plt.ylabel(\"Pressure (MPa)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below are the depths at which there is an overpressure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1300.,  1400.,  1500.,  1600.,  1700.,  1800.,  1900.,  2000.])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "depth[-(len(porosity)):][Pp > hydro[-(len(porosity)):]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the overpressure starts at between 1200 and 1300 meters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem 4 (20 points)\n",
    "\n",
    "Compare the characteristic hydraulic time of escape in a sandstone formation with porosity 0.25 and the one with the same matrix grains, but more compacted, with porosity 0.12.\n",
    "\n",
    " * You can assume that permeability can be estimated using equation $k = C d^2 \\phi^3$ where $C$ is a given constant, $d$ is sandstone grain diameter and $\\phi$ is porosity. When using the formula, assume that constant $C$ and grain diameter $d$ do not change during compaction. Compressibility $\\beta$ of a material is inverse of its bulk modulus $K$. Bulk modulus for brine water and quartz grains are given as $K_f=2.2$ GPa and $K_r=41$ GPa."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**\n",
    "\n",
    "We'll make the comparison through utilizing a ratio of the two characteristic hydraulic times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.1711475531756665"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tau_ratio = ((0.12 * (1. / 2.2e9) + (1. / 41e9)) / \n",
    "             0.12 ** 3.0) / ((0.25 * (1. / 2.2e9) + (1. / 41e9)) / 0.25 ** 3.0); tau_ratio"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Therefore, the more compacted sandstone formation with a porosity of $0.12$ will have a characteristic hydraulic time that is $5.2$ times longer than the sandstone with a porosity value of $0.25$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem 5 (30 points)\n",
    "\n",
    "Assume stresses are given as $S_{hmin}=40$ MPa, $S_{Hmax}$=60 MPa, $S_{v}$=45 MPa and $S_{hmin}$ acts in the East-West direction. For each of the faults below, calculate the normal and shear stress and then determine what kind of fault would it be, if it were to slip.\n",
    " \n",
    " a. Fault with strike north-south, dip 65$^\\circ$ to the east.\n",
    " \n",
    " b. Fault with strike north-south, dip 50$^\\circ$ to the west\n",
    " \n",
    " c. Fault with strike east-west, dip 25$^\\circ$ to the north."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**\n",
    "\n",
    "First we need to order the stresses in terms of the principle stress definitions, i.e. $S_1 = S_{Hmax}$, $S_2 = S_v$, and $S_3 = S{hmin}$.  With this ordering, we can see that we need to rotate by $\\gamma = 90^{\\circ}$ to align $S_v$ down and $S_{hmin}$ with the East-West direction.  First we define a function to calculate the rotation matrix and compute $S_G$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def compute_Sg(S, angles=(0,0,0)):\n",
    "    \n",
    "    alpha, beta, gamma = np.radians(angles)\n",
    "    \n",
    "    Rg = np.array([[np.cos(alpha) * np.cos(beta),  \n",
    "                    np.sin(alpha) * np.cos(beta),  \n",
    "                    -np.sin(beta)],\n",
    "                   [np.cos(alpha) * np.sin(beta) * np.sin(gamma) - np.sin(alpha) * np.cos(gamma), \n",
    "                    np.sin(alpha) * np.sin(beta) * np.sin(gamma) + np.cos(alpha) * np.cos(gamma),  \n",
    "                    np.cos(beta) * np.sin(gamma)],\n",
    "                   [np.cos(alpha) * np.sin(beta) * np.cos(gamma) + np.sin(alpha) * np.sin(gamma), \n",
    "                    np.sin(alpha) * np.sin(beta) * np.cos(gamma) - np.cos(alpha) * np.sin(gamma),  \n",
    "                    np.cos(beta) * np.cos(gamma)]])\n",
    "                  \n",
    "    return np.dot(Rg.T, np.dot(S,Rg))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  6.00000000e+01,   0.00000000e+00,   0.00000000e+00],\n",
       "       [  0.00000000e+00,   4.00000000e+01,   3.06161700e-16],\n",
       "       [  0.00000000e+00,   3.06161700e-16,   4.50000000e+01]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S = np.diag([60, 45, 40])\n",
    "\n",
    "S_G = compute_Sg(S, angles=(0,0,90)); S_G"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll define a function that will compute the unit vectors for us."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_unit_vectors(strike, dip):\n",
    "    \n",
    "    strike = np.radians(strike)\n",
    "    dip = np.radians(dip)\n",
    "    \n",
    "    n = np.array([-np.sin(strike) * np.sin(dip), np.cos(strike) * np.sin(dip), -np.cos(dip) ])\n",
    "    \n",
    "    ns = np.array([ np.cos(strike), np.sin(strike), 0 ])\n",
    "    \n",
    "    nd = np.array([ -np.sin(strike) * np.cos(dip), np.cos(strike) * np.cos(dip), np.sin(dip) ])\n",
    "    \n",
    "    return (n, ns, nd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With this we can easily compute the anwers to the following scenerios."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**a) $strike = 0^{\\circ}, dip = 65^{\\circ}$**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-0.          0.90630779 -0.42261826]\n",
      "[ 1.  0.  0.]\n",
      "[-0.          0.42261826  0.90630779]\n"
     ]
    }
   ],
   "source": [
    "n, ns, nd = compute_unit_vectors(0, 65); \n",
    "print(n)\n",
    "print(ns)\n",
    "print(nd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we compute the normal and shear stresses on the plane"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "40.893030975783645"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sigma_n = np.dot(np.dot(S_G, n), n); sigma_n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tau_s = np.dot(np.dot(S_G, n), ns); tau_s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1.9151111077974452"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tau_d = np.dot(np.dot(S_G, n), nd); tau_d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, there is only a shear stress in direction of the dip.  This is the shear stress on the footwall as defined, therefore the hanging wall will have a shear stress in the opposite direction (opposing motion) such that the hanging wall will move in the direction of the shear stress on the footwall.  Therefore, the motion of the hanging wall will be *up* the footwall.  **This is a reverse fault**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**b)** $strike = 180^{\\circ}, dip = 50^{\\circ}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ -9.38133875e-17  -7.66044443e-01  -6.42787610e-01]\n",
      "[ -1.00000000e+00   1.22464680e-16   0.00000000e+00]\n",
      "[ -7.87187789e-17  -6.42787610e-01   7.66044443e-01]\n"
     ]
    }
   ],
   "source": [
    "n, ns, nd = compute_unit_vectors(180, 50); \n",
    "print(n)\n",
    "print(ns)\n",
    "print(nd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "42.06587955583268"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sigma_n = np.dot(np.dot(S_G, n), n); sigma_n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.8762677505405461e-15"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tau_s = np.dot(np.dot(S_G, n), ns); tau_s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-2.4620193825305208"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tau_d = np.dot(np.dot(S_G, n), nd); tau_d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Following the same logic as above, **this is a reverse fault**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**c)** $strike = 270^{\\circ}, dip = 25^{\\circ}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[  4.22618262e-01  -7.76337152e-17  -9.06307787e-01]\n",
      "[ -1.83697020e-16  -1.00000000e+00   0.00000000e+00]\n",
      "[  9.06307787e-01  -1.66486040e-16   4.22618262e-01]\n"
     ]
    }
   ],
   "source": [
    "n, ns, nd = compute_unit_vectors(270, 25); \n",
    "print(n)\n",
    "print(ns)\n",
    "print(nd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "47.679092927350951"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sigma_n = np.dot(np.dot(S_G, n), n); sigma_n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1.2751975718967087e-15"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tau_s = np.dot(np.dot(S_G, n), ns); tau_s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.7453333233923374"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tau_d = np.dot(np.dot(S_G, n), nd); tau_d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Following the same logic as explained in a), **this is a normal fault**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem 6 (20 points)\n",
    "\n",
    "A major fault in your reservoir has strike N35$^\\circ$W and dip 60$^\\circ$ from horizontal. $S_{hmin}$ in this area acts along 020$^\\circ$, and in situ stresses are $S_{hmin}=70$ MPa, $S_{Hmax}=80$ MPa, $S_v=65$ MPa.   Assume the pore pressure $P_p =30$ MPa. Is the fault likely to slip if the respective frictional coefficient is $\\mu=0.7$?\n",
    "\n",
    "**Solution**\n",
    "\n",
    "First we'll define some functions that we will use in this problem and the next."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def compute_stress_components_on_fault(S, strike, dip, angles=(0,0,0)):\n",
    "\n",
    "    S_G = compute_Sg(S, angles)\n",
    "    \n",
    "    n, ns, nd = compute_unit_vectors(strike, dip)\n",
    "    \n",
    "    sigma_n = np.dot(np.dot(S_G, n), n)\n",
    "    tau_s = np.dot(np.dot(S_G, n), ns)\n",
    "    tau_d = np.dot(np.dot(S_G, n), nd)\n",
    "    \n",
    "    tau_mag = np.sqrt(tau_s ** 2. + tau_d ** 2.)\n",
    "    \n",
    "    return (sigma_n, tau_s, tau_d, tau_mag)\n",
    "\n",
    "\n",
    "def check_fault_slip(S, Pp, strike, dip, angles=(0,0,0), mu=0.6):\n",
    "    \n",
    "    sigma_n, tau_s, tau_d, tau_mag = compute_stress_components_on_fault(S, strike, dip, angles)\n",
    "    \n",
    "    return (tau_mag / (sigma_n - Pp) > mu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can check to see if our fault is slipping under the given conditions from the funcions above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S = np.diag([80, 70, 65])\n",
    "\n",
    "check_fault_slip(S, 30, -35, 60, angles=(-70,0,0), mu=0.7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Problem 7 (40 points)\n",
    "\n",
    "Assuming the same reservior conditions as Problem 6, write a computer code that can compute the resolved shear stress magnitude and normal effective stress on a fault characterized by the strike direction and dip angle.  Now randomly generate 10000 faults (use a statistical distribution to choose random pairs of strike and dip) and plot the results on a Mohr diagram.  Where do near vertical strike-slip faults appear on this diagram? Where do the oblique-slip faults fall?\n",
    "\n",
    "**Solution**\n",
    "\n",
    "The code below creates the random faults.  We will assume that a near-vertical strike-slip fault is one in which the dip angle is greater than $85^\\circ$.  Oblique-slip faults are those that have near-equal magnitudes of shear and dip components of shear stresses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe0AAADfCAYAAADbcYvpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XlclXXa+PHP94Bs7iDuKdaIipJa\njg4u6BO9UCctq6cxB55m+k0DYtNT2ShN02g100Y9T2WlSOuYTjVPmU1NmiOTMY6UWeGKSwuWC4qa\nuKCynO/vj/scOMDZgLNzvV8vXnjuc3Pfl8DhOt/t+iqtNUIIIYQIfCZ/ByCEEEII90jSFkIIIYKE\nJG0hhBAiSEjSFkIIIYKEJG0hhBAiSEjSFkIIIYKEJG0hhBAiSEjSFkIIIYKEJG0hhBAiSEjSFkII\nIYJEuL8DsNWjRw+dkJDg7zCEECJofP7558e11vFtvEbP8PDwF4ERSGPOn8zAztra2tuuvPLKY/ZO\nCKiknZCQwNatW/0dhhBCBA2l1IG2XiM8PPzF3r17D4uPj//BZDLJhhR+YjabVUVFRVJ5efmLwLX2\nzpF3VEIIIUbEx8efloTtXyaTScfHx1di9HjYP8eH8QghhAhMJknYgcHyc3CYmyVpCyGECAhff/11\nh7S0tMsGDhw44pJLLhlx6623XnLhwgW1efPm6DfffLOr9bz58+f3XbRoUS9/xuovATWmLUQgKC4u\nZuPGjWzcuJFNmzYRFRXFiBEjiI2NpXfv3txyyy0ArFixAoDRo0dz4sQJpkyZQkpKij9DFyJomc1m\nZs2a9aPbbrvt2J133vl1bW0tP//5zwfeeeed/YYPH35+69atHWfPnl3piXvV1tYSHh6c6S84oxbC\nheLiYlasWEF5eTk7d+7k2LFj/OQnPyEmJoaSkhJiYmKYMWMGp0+f5pNPPuHQoUNERkZy9uxZTp06\n1ehaVVVVFBUV1T/Oz8+3e0+lFP369aNLly5ERERw8eJFIiMjOXXqFBcuXAAgNjaWGTNm0K1bN0ny\nIqht2LChY2FhYee0tLQzV1999bm2Xu+9997rHBkZab7zzjtPAISHh5Ofn/99QkLC5eHh4VprzdCh\nQzvdc889RwBKS0ujx44dO+Tw4cMRc+fOPXr//fcfA1i6dGnssmXLetXU1Kgrrrji3IoVKw6Eh4cT\nExMzOisr6+g///nPLk888cTBqVOnnm1rzP4gSVsEreLiYvLy8igqKqKqqooePXrQt29fjh49yoED\nzSfUrl+/vtHj3bt3ezQerTUHDx50ek55eXmj+6ampvLYY4+RkpJS38KXZC4C3YYNGzrOmDEjsaam\nxvTUU0+Z33///X1tTdw7duyIHjlyZJXtsdjYWHO/fv2qMzIyju/bty9qxYoV3wHMnz8/+quvvora\nvHnz3lOnToUNGzZsxIIFCyp27doV+dZbb8Vu3bp1T2RkpM7MzByQn58f95vf/ObE+fPnTSNGjDj/\n9NNPH25LnP4mSVsEldzcXF555RVOnz7NxYsXGz138OBBl0kz0BQVFZGamsrs2bNZtWpV/fGFCxdy\n2WWX8fbbb3PjjTeSlZXlxyiFaKywsLBzTU2NyWw2U1tbayosLOzc1qSttUYp1WwynOV4s/PT09NP\nRUdH6+jo6NrY2NiagwcPhq9bt67zzp07Y0aOHDkM4MKFC6aePXvWAoSFhfHLX/7yh7bEGAgkaYuA\nlJmZyTvvvENYWBgdOnSguroas9lMVVWV6y8OMrW1tY0SNkBeXl79v9evX8/8+fO58sorSUpKkjF0\n4XdpaWlnnnrqKXNtba0pPDzcnJaWdqat10xOTj7/7rvvdrc9dvLkSVN5eXlEWFhYs2QeGRlZfyws\nLIza2lqltVY33XTTieeff/5Q0/MjIiLMwTqObUtmjwu/Ky4u5vrrr6dPnz5069aN7t27s2rVKqqq\nqjhz5gwnT57k7NmzIZmw3XXu3DmKiorIz88nOzub++67j/Hjx3PJJZeQk5NDcXGxv0MU7cjVV199\n7v3339+3YMGCQ57oGge49tprz1y4cMH03HPPxYHxZnbevHmX3HTTTcd79+5dc/bsWZf5atq0aaff\nf//97ocOHQoHOHr0aNi+ffsi2hpbIJGkLfxq6tSpjB8/njVr1lBeXk5lZWWziWDCsYMHD5Kfn8/4\n8ePp27cvffr0Yfjw4RQUFPg7NBHirr766nOPPvpouScSNoDJZGLNmjVfrV69uvvAgQNHDBo0aERk\nZKR5yZIlh6ZPn35m37590UOHDk164YUXuju6xpVXXnnh/vvvP5SWlpaYmJiYdNVVVyV+//33HTwR\nX6BQWgfOevoxY8ZoKWMaugoKCnjppZf44YcfOHjwIOfPn/d3SI3079+f06dPc/r0aX+H4hF9+vRh\n3LhxLFy4ULrRQ5hS6nOt9Zi2XGPbtm1lI0eOPO6pmETbbNu2rcfIkSMT7D0X/B38IqDZrnluOns7\n0Bw8eJD09HQKCwupq6vzdzhtduTIEdasWcOaNWtISEhg1KhRksCFCHKStIXHWZdiffnll3aXXgWy\ntr6x6Ny5M2fOtHlOjseVlZVRVlYmCVyIICdj2sJjCgoKSEhIqB+jDraE7QmBmLCbsibv8ePH06dP\nHxn/FiKISNIWbVJcXMyjjz7KuHHjyM7ObpeJui06derk1/uXl5eTnZ1NYmKizEAXIghI0hatUlxc\nzOTJk5kwYQL33XcfW7Zs8XdIbhk8eLC/Q2jEXtEIf9i/fz/jx49n9OjRkryFCGAypi3cZq3n/ckn\nn1BSUuLvcFpl//79/g6hkXPnXK+WMZlMmM1mH0QDJSUljB8/npEjRzJo0KD6DVJk7FuIwCAtbeFS\nbm4uPXv2ZPz48eTn5wdtwg5E7iRjpRS9e/f2QTQNtm3bxpo1a+rXgHfp0kXGvoXXOdqac8mSJXG3\n3HLLAHtfExMTMxqgrKysw7Rp0y71bcQG261C77rrrr5r1qzp7K17SdIWTiUlJZGXl0dFRYW/Qwlo\nCQkJXrt2XV2d39eOnzlzhuzsbK6//nrpPhfkf5Yf2/d/+iabHjRd2fd/+ibnf5Yf29ZrWrfmvPba\na08dOHBg57fffrvz3LlzpjvvvLOfO1+fkJBQs27dum/aGkdbPf3004dnzZrltRmpkrSFXQUFBcTF\nxVFaWurvUHyqc+eWv0E2mUxMmzYNk8kzL6cOHTowcuTIRscCpYTrmjVrmDBhApMnT5bk3U7lf5Yf\ne/f6uwceOXskQqM5cvZIxN3r7x7Y1sTtaGvON998s0dVVZXp0KFDHSZNmjQ4ISFhxD333NOn6dfv\n3bs3YvDgwcMBzp49q2bMmHFpYmJi0jXXXHPp5ZdfPrSoqCgGGlrmAK+88kr3G2+8MQHg8OHD4VOn\nTr1sxIgRw0aMGDFs/fr1He3FOW/evH6XXXbZ8MTExKSsrKz+TZ+/8cYbE1555ZXuAP369UvOycnp\nl5ycPCw5OXnYzp07I9vyPQIZ0xY2CgoKePrppzl16hRHjhzxdzgt4qn10U6vkQykAV0BM8ZbXjOY\nTWb+0uEvmH9rhmigEigEdhhfppTCZDK5VbBl8ODBpKWlAUYXdSDSWlNUVMT48eMZOHAg9913n+xC\n1o48VPRQvwu1Fxq9Q71Qe8H0UNFD/eb+eO7J1l7X0dacffr0qa6trVXbt2/vuGPHjl2dOnUyjx49\nOum6666rTE1Ntftu9sknn+wZHR1t3rdv3+5PP/00esKECUmu7p+dnX3J/Pnzj06dOvXs/v37I6ZO\nnTr4m2++2WV7ztGjR8M++OCD7t98881Ok8nE8ePHw1xdt0uXLnU7duwofe655+LuuOOOSz766KOv\nXH2NM5K0BWDsqtV0p6lg0pKErZSiUfle22RsTbjQPEEDWCd7hzX+fLr2NMRYjnUDZlr+vcNIcjNn\nzuTw4cONZ9knA9No+Loq2L9uP/v/vd+492KavQGwlZqayubNm6mtra0/1rt3b44dO+Z0rDwiIoLq\n6mqHz7fEgQMHyM7OJjc3l6ysLB5//HGPXFcErvKz5XY34HB03F2utuacOHHi6d69e9cBXHPNNT9s\n3Lixk6OkvWnTpk7//d//fQxg3Lhx5xMTE112Vf373//usn///mjr47Nnz4b98MMPpu7du9e/mGJj\nY+siIyPNN99888Brrrmmcvbs2ZWurvuLX/ziJMCvf/3rk/fff/8lrs53RbrH27GCggIGDRpEeHh4\nUCfsFkkGfac2EuJdQCZwA0aiVZbPNzQ5Fmb53JLVWREYidfigw8+YMScEbAQ496LLffoaHPtjpZj\ns5rEM9OI20opRXR0NI899hjz589vtGzs+HHX5aM9lbBtnTp1iry8PGJiYsjNzfX49UXg6N2pt91f\nIEfH3ZWcnHy+pKSkUZe07dacTZdHulou6eh52+Pnz5+vf6C1ZuvWraV79uzZvWfPnt3Hjh3b3r17\nd/PEiRMHDx06NGn27NkDO3ToQElJSemNN954as2aNd2mTJnicg2p7bCZvTclLSVJux0qLi6mZ8+e\nZGdnU1ZW5tE624Gy7phkYAENCXIBMB0jAdomxMtonoxbmqAd6YrxxiAZqq+u5uXKl41WtcLxPaxv\nEmxFYLTI77L8XxbB+dzzTFw/kbzoPONNiCWp19bWNmpl+/rncf78efLy8sjMzPTpfYXvLEpddCgq\nPKpRV05UeJR5UeqiZntYt4SzrTljYmLMmzZt6nL06NGws2fPqg8++KDb5MmTzzq61sSJE8+uXLky\nFuCzzz6L2rdvn7U/i7i4uJovvvgiqq6uDtv9uydOnHj68ccf72l9vHnz5miATZs27d+zZ8/uN998\n80BlZaXp5MmTYbNnz67Mz8//vrS0NAYXVqxYEQvw0ksvdR89enSbd0ST7vF2Zty4cV4thOL3XeOm\nAz+meVLsCIzFfoL2FtuWe1tZkz2gMb7HZsyNW+MjgUHUj7WzFfRa//w8Vq1axe233y7ru0OQddz6\noaKH+pWfLY/o3al39aLURYfaMp4NDVtzZmVlDXziiSf6mM1mrrrqqsolS5YceuGFF2LHjBlzdvbs\n2YPKysqibrzxxhOOusYBfvvb3x67+eabByUmJiYNHz68Kjk5uT5ZPvjgg4euu+66H/Xp06dm6NCh\n58+dO2cCKCgo+P62224bkJiYmFRXV6fGjRt3Zvz48d/ZXvfUqVNhM2bM+NHFixcVwJ/+9KfvXf2/\nLl68qC6//PKhZrNZvfHGG22e3S5bc7YjSUlJoT0bPBP7Lef2QtP4/66Br4GV/gmnU6dO1NTUMHTo\nUJYtWyYJ3Etka07Xxo4dO+TJJ5/83lmi94Z+/folb926tbRPnz61rs9uIFtztnO5ubn87//+b6MJ\nS8HE5cSxvcAIGrVG2yV7vQiX0TAe3nRinZNJbp5w9qzRe7lt2zYmTpzIpk2bJHEL0UaStENYQUEB\nCxYs8HthDoCBAwdy8eJFysvLGx1XStGrV69mx+tNB/1jbSQgDdQCHWhIUN2w3+0tDIqG7vmmM99t\nJ915OYGbzWbuuusunn76aUncwue2bNmy1x/3PXTokMdfUTIRLQRZN/PIzs4OiIQNxtIge4lZa+00\nYTMW47dUWT5H4Ntx6VDgbGKd9bkms9QdTWBrSznVLVu2MGnSJC655BKZYS5EK0nSDjEFBQWMHz+e\noqIif4fiVM+elkmayTTMiv4DDUuxkoExSEL2pQiM5WbJxqSgsLDmdSOOHTvWplvU1dVx8OBB8vLy\niIyMlOQtRAtJ93gICeQCKY3GpZPhWNoxY3wVHHfbCt8LA26AuhvqGuYLDKF+/oC50OyxLvTq6mry\n8vIApCiLEG6SlnYIKC4uplu3bgGbsMFmKdh0GhcucdVtK3zPtst8LM0Lz9xHo2IvbZWXlyf7eAvh\nJq8mbaVUN6XUW0qpPUqpUqWUzEDxsMzMTMaPH09lpctqej7TaDzUtvv7d8iksWBjb/6ATTe6yy93\ns7iLdR9v2f6z/VJKXfnrX/+6fgOORYsW9Zo/f35ff8bkyL333ttocsfo0aOHtuY6tlt6usvbLe1n\ngHVa66EYpR9CeJGwbxUUFNCnTx+Ptq7DwsI8UkHLthu8UQWySCRhh4owGpVptVJKkZqaytixY0lI\nSGhxsZ3s7Gw6duwoY92BLj8/lr59kzGZrqRv32Ty2741Z0REhP7ggw+6HzlyxKPDtjU1NR67ltls\npq6ujiVLljTaZezLL7/c47GbuOC1pK2U6gKkAi8BaK2rtdanvHW/9iQ3N5fs7GzHs65bQSlFXV1d\ni/7IxsTE0K1bN8cnpGG0yoTX5IzJwaT8NMplmZOglKrvUdGLNCVTSviy9kvKyspaddmqqiophRrI\n8vNjufvugRw5EoHWcORIBHffPbCtiTssLEzfcsstFY888kizlqejbTM/+uijmNGjRw8dNmxY0ujR\no4du27YtEmDJkiVx06dPv/Sqq6760aRJkxJtr5WTk9Pvsccei7c+nj9/ft/Fixf3AvjDH/7Qa8SI\nEcMSExOT7r777r5gbPl56aWXDs/MzBwwfPjwpNmzZydcvHjRNHTo0KRrr712EDTe7vP+++/vlZiY\nmDRkyJCkefPm9QP4n//5nx4jRowYNmTIkKSpU6dedubMmVa/aL35ar8UqABeUUp9qZR6USlld39S\n4Z7i4mJGjx5dP3nHk1pTGa+qqsr+krLpGDPBuzZ/SnjWX3f9lf9I+A//3FwDi4xEbTtP4bTpNDU/\nrWnzuPeqVaukxR2IHnqoHxcab83JhQsmHnqoX1svvWDBgmOrV6+OPXHiRKOlC9ZtM3fu3Fn6zjvv\nfD137twEgJEjR17YsmXLntLS0t2LFy8+tHDhwvru9S+++KLT66+//u0nn3yyz/ZamZmZJ99+++36\nNxjvvvtu98zMzB9Wr17d5auvvoravn17aWlp6e6SkpKYtWvXdgIoKyuLuvXWW0+Ulpbufuutt8oi\nIyPNe/bs2f23v/3tW9tr//Wvf+3y97//vfvnn3++Z+/evbsXL15cDpCRkfHDzp07S/fu3bt7yJAh\n55csWdKjtd8jb84eDweuAO7QWn+qlHoGuBfjz3k9pVQWkAUwYMAAL4YT3AoKCsjOzvZ3GM3Ub05h\nW6UMpBvcR06cP0Hht4WuT/QGZ2/5w2ko3FIFrKNVs87z8vJYv349S5culaIsgaLcwRacjo63QGxs\nrPmmm2468dhjj/WMjo6u35TE0baZls07BpWVlUUppXRNTU39X55Jkyad7tWrV7PdkCZMmHD+xIkT\n4WVlZR2OHDkS3rVr17rBgwdXP/HEEz2Lioq6JCUlJQFUVVWZ9uzZE3XppZdW9+nTpzotLc3lZh//\n+Mc/umRmZh7v3LmzGcB6/88//zx60aJF/c6cORN27ty5sMmTJ7d6EpI3k/ZB4KDW+lPL47cwknYj\nWusCoACM2uNejCdo5ebmeqV13Sa2Sdpa81oStbBl/X2wbjk6jVYl75KSEiZOnMiyZcvIysryaIii\nFXr3rubIkeYJunfbtua0+t3vfnf0iiuuSLr55pvra6Fbt83s1KlToxxx2223DZg8efKZf/zjH1/v\n3bs34qqrrhpifS4mJsbhpvIzZ878YeXKld3Ly8s73HjjjSet97jrrruOLFiwoFEN9r1790Y4u5Yt\n697fTWVlZQ166623vkpJSTm/ZMmSuI8//rizO9ezx2vd41rrcuB7pZT1m5gG7PbW/UJVZmZmYCXs\nZIy3XrbLtqwVy4RwxHa/8AW0uOvcbDYzb948WRYWCBYtOkRU4605iYoys6htW3Na9erVq27mzJk/\n/OUvf6nvQna0bebp06fD+vfvXw2wfPlyt7uc/+u//uvk22+/Hfv+++93z8zM/AFg+vTpp1977bUe\nlZWVJoBvv/22w6FDh+w2bMPDw7V1py9b06ZNO/3aa6/1sI5ZHz16NAyMVvuAAQNqLl68qN544402\njf17ewbLHcAqpdR2YBTwiJfvF1IKCgoCa+11MnAdEIUkaWDlDSuJ6eByO11hqw3Ju66ujgceeEAS\nt7/NnXuSp546QJ8+1SgFffpU89RTB5jbtq05bf3+978vP3XqVH3CLCgo+P6LL77omJiYmHTZZZcN\nf+655+IBcnNzyx944IH+V1xxxdC6umY94Q6NGTPmwrlz50y9evWqHjhwYA3ADTfccPqmm246+eMf\n/3hoYmJi0vXXX3/ZqVOnmpcFBDIyMiqGDRtWPxHN6j//8z9PT58+/dSoUaOGDR06NOmPf/xjb4B7\n77338NixY4dNmjQpcfDgwRda8S2pJ1tzBqDi4mJ+9rOfcfDgQX+HYpiOUVJUWtT1wlQYf77+zwBk\nrpZZzq1WDbxHi7rMw8PDef7556Wr3EK25gw9zrbmlIpoAaa4uJjx48cHRsKeDizCKIgShiRsG3W6\njszVmZKw2yoCuBrnSwebqK2tJTs7m8mTJ0urW7Q7krQDSEFBAVOnTvXJvVwWUWm6w5bwP235CDVd\n4NQpo4RDbGys2wV+ioqKmDBhglRRE+2KbBgSIMaNG8eWLVt8dj+7wyKybCuwWfe8vg77r1zrLH5H\nrD/yQPy5LgYq4WThSWjByKjWmuzsbL7++mvZdES0C9LSDgBTp071acK2y92NPIT32WtN12Ik7B3A\nu8A5Glre2vJ4C3AK+y3yamA1xprpQGO7QckNNN6e1U15eXnS4m4bs9lslld9ALD8HBwuMZOWtp8V\nFxezfv16/9xcWtb+UQdcAKwTz22/79XAl8AIm+ebFCcx7TKhdxpZuVmPyVrLZ9ufrbWFbp3s5ail\nDq5b695mvbc1gc8A3setiWr3338/ycnJUoSldXZWVFQkxcfHV5pMplAchAkKZrNZVVRUdAV2OjpH\nkrYf5ebm8tRTT/nl3upyhZ6hpTa4N9hLfGbLsUroW9qXwx8eNo47Sq5rcai+Cp0zO4wPpRRKqYav\nsSa/a2i+gUs1UIOxJCsQWDeZmWV57CJxV1RUMH78eBYuXChd5S1UW1t7W3l5+Yvl5eUjkB5YfzID\nO2tra29zdIIs+fITX49hN7OAwPnjHCo0RuIrAYYAXSE2LJY0lcbqB1ZTV1eHUoro6Giqqhr3Uyck\nJLR6gw1b0dHRnD9/vtExpZTrOQzWNwxg7MwWaG/mzgFPuH/68uXL282SME8s+RLBQ1rafuCvMeyR\nI0eyU+2kbnpd4P1RDjbWceNqjNZg0y5oS0t5wKgBFB83liVZk2fThD148GA6dmzbOyjrtZsmbKd2\n4Lj16mjYxF/d5zHUT1Zr9H124De/+Y1MThMhSVraPpaZmemTKmcJCQlER0dTWtqwhXn8L+OpSKjw\n+r1DXisKgoD9VrCVyWRq1u09a9YsysrKKCkpaV2cFv3792/9un/bwjpm4FugP/7dG70WYzKeG9//\nsWPH8umnn7o+MYhJS7t9kbELH8rNzfVZWdJz585xxRVXNByYDhUDJWE7ZTsb28Hzkecj7SbsHj1c\nlz121gpumrCta5W3b99u93x31zIDHDrUhpLQayEyL5JZ22ZhetgEK4HHgK/x35rxcIzNR9ywZcsW\n2d5ThBRJ2j7i6526KioqGt4gWAulyOxw5yqBBzGWRtU2fkqZFTOqZzB973RiDzWu9z927Nj64iBW\nJpOpRYm1KZPJxLvvvmt30llYWFiLrt3W3rSamhoSExPp0KFDw8GVGN+npkvPNMbMeG8n9BjcXhYm\nS8FEKJHucR+YOnWq/5Z1ScJ2nwbTH433sebhZsdLptxkr8vbXQ4nj9kxB2MnngHAd8B9wHhgLkb1\nWaszGHMPrT3d54BOwFmb42Dk2xOWf8dZrvl74MthwxoNtzg1F+iF93/v3ByqiI+P59Zbbw3JMW7p\nHm9fJGl7ma/GsJux/tEESdjuOgUdnu/A8OHD2zyO7A3PAjnY7x5rOles6bG2sq5YAzgO3Am87uqL\nMoHLvBBMU6eAp907NRSXg0nSbl+ke9yLiouL/ZuwpbKZY/YqhhUam1Fs27bNI7cwmVr38pqDMd+r\nDqOX3mz5uJ2GfVuaftjyxo/dWoJeAfHAKpu4Ki0xW82aNYuFCxeiVqmG4QZHldo8oSv13eSuhg2W\nLFkim4yIoCYtbS/q3LkzZ8+e9f2NFyPJ2hlr1THLWmp3ur/Hjh3LsWPH3F5LbTKZWLZsGU8++ST7\n9+93O7RnMZJzsP34bP+KmIF84Df2TrS2vj39H7StMufi52kymdi0aVPIVE6Tlnb7IknbS3xePEVK\nkjpnU5HM3h90k8mE1trhOHJYWBhaa7fHqFNTU/n444/p2bMnFRX2Z+3PwUhunZscD4UfXdPvoqZh\nvP11225z8M5/2MWyMOvPJxRI0m5fpLiKFyQkJHDgwAHf3TCZwKxiFUgURletHdZubGdvYOvq6urP\njYuLc5iIrT799FNyc3Prl3nZThbTNIxLhUKCtsdel30CRrf6qpXGsfXANNs3m9bCLZ74poRjlGp1\nkLS/+eYbD9xECN+TMW0Pi4uL823CBuOPniRs5yqbHzKZTKSnpzdqQSulnI5Fh4eH86c//YmwsDCH\n5wBUV1eTl5fHo2fPUoeRrBIwXnC249Ltje04fDpwYQfMeRrjDdVDeHYXskiMcr12loVVVVVx/fXX\ny/i2CDqStD1o6tSpnDzZgs2A2yoZY61qV1cntnOWSWZNaa1Zv359sxb2nDlzGq9JtmFtcTvrJp8D\nnNK6fvKYdRKXz0W48U7OxZsPb7LuB7IKY0j6WTB2M/PUiJ3CWMt2A8ZYuo2TJ0+yZs0aJkyYIIlb\nBBVJ2h5SUFDg27XY1i5x6/7XokEdDUU/TuFwHa+97nCtNW+++SZ1dXV2W9N1dXXk5uY2T/Q0zPhe\nBXTBR63pnBzQ2v7HxYuOn7N+1NbaP75yJcTFNb5XRAS0cka8MwrjD9HtgHkHlDyJZ2eZK4wx9Mzm\nT2mtycnJ8eDNhPAuSdoekJmZSXZ2dou/rn///q2/qXSJ26eBNRg7Qj2IsX63hUVRamtrMZvN9a3q\npppWP1tH4+5vjydqa6IcONBIprbJdelST9/NkJEBx483fxNQV+c4qbeR9U3O5efA/CCULMGzre7L\nMIoNNbF9+3apmiaChiTtNipUnd4WAAAgAElEQVQoKGj1WuwjR460/sbSJW5fJS1O0q1h27JOxwuJ\nWqmGVrQ1UZaVGck0ENhL6itXQht3KwOb5H3SSN7mB6DyEZhjvwx7yy48lmZj3FprsrOzJXGLoCBJ\nu40WLFjQ6q911JJzyY16y+2SZeza1SQxq4SEBIYNG9bi28wBVuCFlrVtS9ps9l4r2lsyMuDsWSP+\ntLQ2X8520lqXali1Gtb92QMXdbDZyEsvvdTGiwvhfZK022DcuHGcPn3a9zeehoxjN1UHvAdqp2Lm\nzJmkpqa6/JKf/exnTJ482a3NN5qOWbd5rWR4uNGatk3UgdSSbqsNGxq3wAcObPMlFZD+rdHyXvdq\nGy4Ug91u8i1btkhrWwQ8Ka7SSgUFBa0ax/YIqXjWmKWQhtqpiIiIQGtNXV2d02Ioo0aNYu/evVRX\nVzvs8bDuQ+3RKmUmE2RnB18r2tPmzYP8fCOpt4IGtsfBqDtaeX8NbAHWNn8qIyODlStXtvLCvifF\nVdoXaWm3QnFxsf8StmjMTH3lq+7duzN9+nTq6uqoq6tzmLAjIiJISEjgwoULTocoJh08yDHamLA7\ndWo8eayuThI2GN8Ds9n4nvTt2+IvV8DlJ9ow3m0d315Ms1b3qlWrmDp1aotjEsIXJGm3Qnp6un8D\n8GQBimBWC7xD/cSzH374gbVr1zotjpKQkMBPf/pT/v73vzutgGadER5PGxJ2Tg6cORM6Xd7ecuiQ\n8camhfuPNx3vfvb9Ft7XeoGxNEvc69evJzc3t4UXFML7JGm30NSpU/2zCYitnXhnt6RgYV1/3aS2\ntNaa2tpaRo8e7fBLy8rKWLNmDTU1NQ7PeZY2zgiPizOSkLSo3ZeRYbS8W7mUTAG3b21F4rZ+8Via\nrePOy8uTMW4RcCRpt5BPC6jYkwyMpv2OaWuMrR5t1l9bJ5KZTCYiIiKYMmVKiy87BzhGwxaYLf72\nhoU1dIMfPy6t69ayLiVrRfK2Jm7zA3Ds8RZ2mTsowJKTkyMV00RAkaTdAgkJCf4OQYqqVNFsHXZ4\neDjp6elkZWXxi1/8gtOnT7s1I9zqWRp3hbc4YefkGJXFJFF7jm3ybsHa7/o9v8/Dq6tbmbhtmM1m\n8vLyWnARIbxLZo+7KTc31ycvXpPJ5Hz7x/Y8c9zaynayT7LWGqWU21totmlm+MCB8PDDkqx9YdUq\nY9b9uXMt+jINHI+GO6fD65e7+QV2fsc2b94csPtvy+zx9sXrLW2lVJhS6kulVGtGmwJCcXGxTxJ2\n586dXScb93JR6LEu0XFS7cxsNru15/UcjMJpreoKDw8PzXXVgc5auGXlSvc2QrGwtrpfedfNVrfC\n2GDkXhoVMZo3b17L4hXCS3zRPX4nUOqD+3iNr2aLnzlzxvVJ7W1AQwMXMVo/dtbUttQc4DVauaFH\nWBi8+qokan/KyDBqoLew2zyyDp5x9/dHAVHALOoTd0lJicwmFwHBqylAKdUfYyv6F715H2/KzMz0\n/2xxW3b2hQ516jFFeGk4YWFhhIWF0aVLl1ZdZw6wEmM/6xZLS5Nx60BibXm3oFxqj/MtHOMOo1HJ\n0yeeeEImpQm/83a77WlgIUHcqdvazUC8ppD2tdyryljKNWLECCZMmEBdXZ3bpWMjbLpRrZPNWvQL\nHxXV0BW+YUNLvlL4yoYNxkRANyYeKoz13OYHoOQ5N68fQ31rW2vNvffe28pAhfAMl3/DlFKfKaVe\nUkrdpZS6SikV786FlVIzgGNa689dnJellNqqlNpaUVHhZti+EefhrQc9wgc7WAWMWowqJxjdk//6\n178aPd2zZ0+nm4NUV1cDrZxslpYG589LyzoYWKurubEvdv0OYseh5kE3Wt4KY8WGRVFRkazdFn7l\nTsPjOuD/MBYazQXKlFIH3Pi6CcC1Sqky4A3gKqVUs4K+WusCrfUYrfWY+Hi33g/4RGZmJidPnvR3\nGPa1hy5yjd3iKbaGDh3KzJkznV6mxQm7Y0ejdS0t6+CzdKnb67sVEK6NlrfLVneTbXDnzp0r3eTC\nb1wmba31Ya31Oq11ntb6Z8AY3Bij1lr/TmvdX2udANwM/FNrneniywJGwHWL2yrE2IbSVqh1mbux\nL3ZsbCz79u1z+PwcWpiwc3KMcVJpXQcv2/XdMTEuT7e2ul1u+bkYuAtIlm5y4V/udI8PsH2stS4F\nhnstogCQlJTk83t269bN5Tn1BUN2AO9hlPK0lvQMJZZ9sZ1RSrFmzRp2795t93nrntctSthSdjR0\nZGRAQYFbW4Jat/x0WALV2qfeDaPfMdnoJpfWtvAHl8VVlFLFwCUY2wnvAC4AV2mtR3k6mEAoruKv\nLTeVUk43sHDpD7RyWnQA0RgVz9YZ22yaTKZGu3DZbgTibC32HOAVINKde3bsCMuXS+s6lA0fDg7e\n3NnSwPNj4I4ZLk48BzwBqampfPzxx56IsE2kuEr74k73eIrWuj9wK/APYBfg6tc6aC1evNgv921z\nZbpQWL9tBp4AdhjfjyuvvLLR01dffTXXXnst0dHRDi9RgjFL3K2EnZQk3eHtwa5dbm3/aa1d7nKr\nT0uvu0xKE/7gTvd4nFIqB5gEHAL+orU+6PXI/KCgoIDy8nJ/h9E6oTA5rclv42effdbo8fr161mz\nZg3nHJSyPAZcjptd4t26GX/MRftw6JDbs8u7VMNr77hI3JatPG+//XbpJhc+5U777B2MvRQewWgH\nVSql9ng1Kj+5/fbb/R1C69mbnBbkWtL7sA7ogZsJOy0NfvihlVGJoLV0qVuJGyBMw0vvOnjSZg/u\n2tpaVqxY4akIhXDJnaTdWWv9EHBUaz2ZhiHDkDJu3Dhqa2v9HUbrWSennSN4Z5JXNfxTKeX2Tl1u\n738dFyeFUto767IwN363ouqczCq3Ju5keOGFF6S1LXzGnaR9wfL5olIqWmv9NvBTL8bkF1u2bPF3\nCC2Wnp7OQNvZsTsw+kJWY7S6gyl52xRSAep363KlReuwn3mmdbGJ0JKRAa+9ZtSSd8I6q9zhOm5L\n4ZW6ujpZAiZ8xp2k/aRSKhZ4E3hZKXUH0M+7YflWZmbQLB9vZMOGDRw8aGd6wQ6MwYyq5k8FJDuF\nVNxJ2OtoQcLOyZEJZ6JBRgb8+c/QqZPT01yu47YUXpElYMJX3Jk9/rbW+qTW+n+BDzCWf13n9ch8\npLi42GeFVNzt7rUa6GKNqdlsbrQkqpl1QE2Lbuk/TQqp/PznPycyMpKwsLD6z7bc7hIHYwxb1mCL\npjIy4MwZl5uOWFvcDlkmpcnYtvCFFi0U0lq/prVeqLUOmWm3vuzWaumyrgMH3KkW68QOCPt7WEMR\nlkDtLrfTI9C5c2cKCwv54x//yEcffcQ999xT/1yLKp2lpckYtnBuwwa3dguzuxTMZlLaSy+9JK1t\n4XWhsLq31QoKCigqKvJ3GB5jryVfV1Jn7LX2oO/jcUuTsWyrV14x5jr+7ne/IyUlpb5inHUWpFsJ\nW2qIC3e5SNzWpWCrVtupnKaAH0NNTY2MbQuva9dJ+6WXXvJ3CB7lqCVfn8wDaS23tfxqk7Fsq5qa\nGjZu3EhxcTE5OTmsW7eOObSgcEpamoxhi5bZsMEouOOEAuZtddDixhjbDtY5MiI4tOuk/fnnTncN\nDRn1ybwQo2XrZ3HRcTwS+YjRA+BkU5C4uDhSU1PJz8+nX1ERq3Czhd23r7SwRevs2uWyq9wEvLLa\n8fN/+ctfpJtceE27TdpTp051PokrFO3AaNn6eS33yfMnmTJlChEREQ7P+e1vf8uXX35Zv3Y+HzcT\ndlKSUf1KiNbasMHlcrAI7LS2k41PWms2btzojciEaL9Je/369f4OwT+sa7n92FU+oOsAUlJSuOuu\nu5o9p5QiNTWVWbNm1ZeULQE6u3PhtDQpTSo8IyvL6dMKeGZtkwM3UD+TfJf8HgovaZdJe9y4cf4O\nwa9SU1P9VvY0pkMMD6c9TEFBAU8++WSz500mE5s2bSI1NZWTJ09Sgpv1xOPipEtceI4bJU97nG9y\nwKZK2qpVq6SLXHhFu0zawVj9zFNMJhNJSUlEfxVtlD31cTf5tF7TeOK/niA7O7vZ9ppKKerq6jCb\nzdTW1nLJv/7l/gYgUu1MeNrSpS7Ht799qkk3uaVKGvh2OaloP9pd0s7NzfV3CH6jlCIyMpJbbrmF\nwsJCEs4k+DyG1QdWs828za1zl2rt/tIumSkuvMHJ+LYCEiphZdNlYF2BZCgtLfVFhKKdaXdJ+9ln\nn/V3CC2ilMJkat2PqXfv3phMpvpknZ2dTWFhISkpKQDGFpe+LnUaBsykftKOLdsla3Nwcxxba0nY\nwrv+/Gfo0MHh0yaMfbjrE7cCrofaYQGwVEOEnHaXtM+fbzoQ5V0tLV3alNa6WTeyu8rLyzGbzYSF\nhbFkyRKWLVsGQE5ODqmpqVRUVBiFTXw9iT6C+i5Ee+YAK3CjW9zJ7HMhPCYjA155BZyUFVYYibu+\nq9wEP1zxAwUFBT4JUbQfqqWlNb1pzJgxeuvWrV67flJSUlB1WUVFRXHhwgWHzyul3C6NGh8fT3h4\nOEePHm3+JiAZmAbEWC/cqnBbRuOwSlsl0MWda0i3uPC1hARwUl64Boh4wPJAQ3pxOh9++KFXQ1JK\nfa61HuPVm4iA0a5a2sGUsAGHCXvUqFEMGzaMSZMmud2Sr6io4MiRI/Zb7dZlYL4sdepgyZnb3eKS\nsIU/PPyw0724w7Fpbdcar1UhPKndJO1QmoBWUlJCaWkpRUVFLd6ExCVfrN+uxlhyZseLuNHQl4Qt\n/CUjA+bOdfi0Ah6x/m6Hw1MnnpIucuFR7SZpP/XUU/4OwedMJhNxcXEt+6JCvLcMTGNUY3sPu+VL\nS4BoV9eQmuLC31xs8zrA+sZXQU3/GrKfy5bELTymXSTtgoICamqCZWNpzzGbzS4n3jXrXndSC7zN\nFMagn517zMGNIiqyzaYIFE42Fvmuq80DBUyDxYsXez0k0T60i6R99913+zsEv6mqcr2mq9mSMm92\nkXeFyMhIunRpPNXsVVwk7JwcSdgicOzaZWxM04QZo6XdqOhKDBy98qhPwxOhq10kbXcSV3tld0mZ\nN7vIK6G2tpbTp0/XH3oWcLwK1sJFl6QQPnfokDG/YuBANEbCNlk+Eirh1dWWxK1Aj9E88PYD/oxW\nhIiQT9qhNAHNZ3YAW/B84tbAXprtrjYPF63sqCgPByKEh2RkQFkZqmPHZn9MI4AX11geKHjwS18u\nzxChKuSTdrBVQAsYa4HVwCmMZOuJBK6AIY0PleDGbPEXX/TAzYXwonPn7B6ONsO6P1sehENmZqbv\nYhIhKeSTtq8roLVWly5dAm9N5w7gaYwxbk8VXOna+KFbk89ktrgIUgpI/7ahxOnf/vY3v8Yjgl+4\nvwPwpkDbGk8pRX5+PidOnCAuLo61a9dy+PBhfvWrX5Fl2b+3uLiYFStWsHv3bioqKoiPj+f06dMc\nOHCAqKgoOnXqRGVlJRUVFZ5fo+1IV9enuM1mkltERARUO9kfNDxcJp+J4GAygYNywwqY+znc8VOj\n3n9xcXF9/X8hWiqkk/a8efP8HQKdO3cmLCyMyy+/nMcee6zRi9WaqG2lpKS49YIuLi5m48aN7Nq1\ni/Xr12M2mzlx4oRHY69XCXTzwHUsY9pWf3OWsAFefdUDNxXCB7KzwVLb354wy/tr83AzGzdulKQt\nWi2kk3ZJSYlf7quUonv37jz66KN2E7Mn2EvuBQUF5OTktHqDEYcKMXbmauv+HNYx7bXGPiXpuOga\nl25xESyWLoV9+6DQfqm/OoUxGDkTvor+yqehidASsmPavq5A1KFDB1JTU9m8eXN9q9dbCduRrKws\nNm3axNy5cxk8eLDnLrwD1N89NKht6Wp3mbA7dvTM/YTwlQ0bIC2t2ZxNDeRfaXkQAa9XvO7jwEQo\n8VrSVkpdopT6SClVqpTapZS601v3smfBggU+uU+XLl1Yvnw51dXVfPzxx37v9kpJSWHZsmXs27eP\n5cuXNyti0lp6m/bMDPJKY122S8uXe+BmQvjYhg3s/Vkatcp4udQqeH4M3DGj4ZTzHYJjcqwITN5s\nadcC92ithwE/AW5XSjmu/edhtsU7vCU9PZ3Kykqft6jdlZWVRWVlJZs3b2bWrFlt3tubtva6V8Oc\nt+B2XLSy4+Kka1wEraFvbuBP/7cY0x+gw+LGCdtq3t/9P99GBCevJW2t9RGt9ReWf58BSoF+3rqf\nLW93jf/oRz9i8+bNXt8n11NSUlJ45513+Pe//82sWbOIjIxs3YW20vrWtgbeg2cOurF67JlnWnkT\nIQLDAzc+wJA9Q5iz3ShpWveATWlTBcs+czxpTQhnfDKmrZRKAEYDn/rifo8++qhXrjt48GA2b97M\n/v37/d4N3hrW5H3hwgXS09NbfoG1bbh5JbADerg6LydHWtkiJPwt6WZeeM8oaWotbfrCew01yX/1\n9K/8Gp8ITl5P2kqpTsDbwF1a62Z91kqpLKXUVqXU1oqKCo/c8/vvv/fIdaw6dOhAeno6+/btC8pk\nbc+HH37I8taMG7emvKkGCo2dvFySGuMiRCS++iodm2wu2LHGst+2gpcrXvZLXCK4eTVpK6U6YCTs\nVVrr1fbO0VoXaK3HaK3HxMfHe+S+TWtbt1ZsbCybN2+muro6aLrCWyIrK4vNmzfTr18LRi3WAl/T\n8sS9A/JxYycvIULFd9/ZPVy/37bLXXKEaM6bs8cV8BJQqrX+X2/dpylP1vZ9//33Q6Zl7UhKSgoH\nDx5k4cKFhIe7uWx/JS1P2snQ2dnzERHSyhahZcAAu4dt99uWnb9ES3mzpT0B+C/gKqVUieXjp168\nHwBvvPFGm68RHx/P5s2bQz5h23r88ccpKipqvre2IzWuT6mnYM4gF+e8LF2FIsQ8/DC1UY0rEp3r\nAPelWR4o+NP2P/k+LhHUvDl7fJPWWmmtL9daj7J8fOCt+1m1tWs8PT2dY8eOtauEbZWSksKmTZvc\n27ikhdXRHvlYqp+JdiYjg/AXX+ZoXBRmoKwr/HomvH55wyl1yjNDeaL9CNmKaK2Rnp4ekmPXLZGS\nksKXX37penZ5pfOnmxrQwvOFCAkZGfQ6fp6wB2DQ3Y0TthCtEVJJuy3rs+Pj49t9wrb14YcfsnDh\nQscnFAIu9vuwdSLayZNpaU6eFCI0PPs+1DwI5geMz9btOlftWOXXuERwCakNQ+67775Wf+27777r\nwUhCw+OPPw5AXl5e8yd3WD5PA2Jw2vc9Zzt0cZTgk5Jk+00R8p59H27f2vAyCdfGY4B7ou8hI1mG\nh4R7Qqql3dqtKZcvX94ux7Dd8fjjj5PhaLx5B8aENBclzh4phEh7Q3dxcbBrVxsjFCLwzf28+ctE\nATmfw9HzR/0RkghSIZW0W2P58uUBWzs8UKxcudLx5LSu9g/bcjie7a39v4UIMGEOlkiaNNyyLcy3\nwYig1q6T9sKFCyVhu2np0qWEhdn54+JqgpmzzcHsXU+IEKQc/K4r4KGPZAa5cF+7TdoZGRn1Y7bC\ntZSUFP71r38xduzYxk8U4rTQypwdTn7JPFS5ToiAl5Xl8GVySaXs+iXcFzJJOy4uzu1ze/fuzcqV\nK70YTWhKSUnh6aefJjraZir4DsfnA7z4rpMh74EDPRWaEIFt6VLOOShb+l1XyN+a79t4RNAKmaR9\n8uRJt8998MEHvRhJaEtJSaGwsJCkJJut0Z1UR4t21ph++GGPxSVEoCt5KKdZ4rZWSNNoWfol3BIy\nSdtdAwcOlHHsNkpJSUFrm86+lpQ0tSVV0EQ7MvHepWTNNCqj2auQdufaO/0anwgOIbVO2x1lZWX+\nDiEkDBkyhNLSUuNBjH9jESJYdP1/OQxKXmZ3zOjEeVlNIVxrVy3tVu0fLexauHAhxkZuOJ1Bvn6Q\ng3lqUgVNtENLr1nqsq6BEM6ERNIeN26cy3O6desm3eIelJKSQn6+ZfJMIVDb/Jw522GIZapBo8Sd\nliZV0IQQohVCImlv2bLF5TkffOD1DcbanaysLKPoyg7gXRpl5jnb4ZV3IaHSaFgoMPbMXrlSErYQ\nQrRSSCRtV4YNGyZlSr0kIsKyR+cO4GLD8WfW2ildWl0Nd8pkGyHmbIdvn4K6B4zPc7b7OyIRLNrF\nRLTdu3f7O4SQ9atf/aqhpyOy4XiP8w6+QEqXinZuznZ44T3oaFl1kVBpPBbCHSHf0o6Pj/d3CCEt\nKyuL1NRUf4chRNB4/J+m+oRt1bHG2FhH1moLV0I+acuWm9732GOPGXXJq/wdiRCBr3+l2e7xAZXw\n+8Lf+zgaEWxCPmnLWLb3paSk8Otf/xrWUT8ZzexoWYsp5H/lhHBKDbBfvve7rvBd5Xc+jkYEm5D+\nC9q7d29/h9Bu3HLLLY3qkG9IaL4+WwNkZ/suKCEC0cMPU+WgnGlsdKx/YhJBI+iTdkFBgcPnpMa4\n76SkpNCzZ0/AmGgz+UDjGhIaUGlpsHSpX+ITImBkZHDXDR0dljMVwhnVqIa0n40ZM0Zv3bq1RV8T\nFRXFxYsX7T4XSP+39iA3N5e86DyO5UG8vdnjcXFw/LjP4xIi0JgeNKHt1ApUKMyL7Y95O6KU+lxr\nPcZTsYnAFvQtbUcJOyZGCmL7mnV/clnuJYRzA7oOsHu8c3hnH0cigk3QJ21Hzp075+8QhBDCrofT\nHqaDqUOzIiszvjgty76EUyGbtIX/yKCEEM5lJGfwy10RvPCeUVzFhPG54D349HGpGigck6QtPE42\nMRLCtfvWnbNbZGX++zKMJByTpC08SzK2EG4Z4GBLW0fHhYAQTdrp6en+DqHdUiiORzt4Mi7Op7EI\nEciq+th/PTg6LgSEaNL+8MMP/R1CuzV3zFzunA4Xwxofrwk3wTPP+CcoIQJQpyeeoTYqotGx2qgI\nOj0hrxPhWEgmbeE/S69ZyrbLkrh1JvXFIw51CafDqysgI8Pf4QkRODIyCH/xZS706oUGLvTqRfiL\nL8vrRDglSVt43IvXv8gb7yoGVUK4Uny3rkj+EAlhR/Gll9Ll5ElMQJeTJym+9FJ/hyQCnFeTtlJq\nmlJqr1LqK6XUvd68lwgczz//fH01Oq01zz//vJ8jEiIw5eXlUVNjTCGvqakhLy/PzxGJQOe1pK2U\nCgOeB6YDScAcpVSSt+4nAsfatWudPhZCGA4fPuz0sRBNebOlPRb4Smv9jda6GngDuM6L9xMBYvr0\n6U4fCyEMU6ZMcfpYiKbCvXjtfsD3No8PAuO8eD8RIFauXAkYLezp06fXPxZCNNatWzeUUmitUUrR\nrVs3f4ckApw3k7a9MhvNKlwqpbKALIABA+wX0RfBRxK1EK5NmTKFqKgoqquriYiIkJa2cMmbSfsg\ncInN4/5AswEbrXUBUADG1pxejEcIIQJKSkoKhYWFbNy4kSlTppCSkuLvkESA82bS/gwYrJQaBBwC\nbgZ+7umbWLuVbB8LIUSwSElJkWQt3Oa1pK21rlVK/Qb4EAgDXtZa7/LSvbxxWSGEECKgeLOljdb6\nA+ADb95DCCGEaC+kIpoQQggRJCRpCyGEEEFCkrYQQggRJFQgTeJSSlUAB/wdh0UP4Li/g3Ah0GMM\n9Pgg8GMM9PhAYvSEtsQ3UGsd78lgROAKqKQdSJRSW7XWY/wdhzOBHmOgxweBH2OgxwcSoycEenwi\ncEj3uBBCCBEkJGkLIYQQQUKStmMF/g7ADYEeY6DHB4EfY6DHBxKjJwR6fCJAyJi2EEIIESSkpS2E\nEEIECUnagFKqm1LqLaXUHqVUqVIqxXL8DqXUXqXULqVUXqDFqJQapZT6RClVopTaqpQa68f4hlji\nsH6cVkrdpZSKVUr9Qym13/K5e4DF94Tle7pdKfWOUspvGxo7itHm+d8qpbRSqkegxRcorxUnP+eA\nea1Y4rzb8r3aqZR6XSkVpZQapJT61PJaeVMpFeHPGEWA0lq3+w/gz8Btln9HAN2A/wA2AJGW4z0D\nMMb1wHTLsZ8CG/39vbTEEgaUAwOBPOBey/F7gccDLL50INxy/PFAiK9pjJbHl2BsvnMA6BFI8QXa\na8VBjAHzWgH6Ad8C0ZbHfwV+afl8s+VYPpDj7++hfATeR7tvaSulugCpwEsAWutqrfUpIAd4TGt9\n0XL8WADGqIEultO6Yme/cj9JA77WWh8ArsN4w4Hl8yy/RdWgPj6t9Xqtda3l+CcY+74HAtvvIcBT\nwEKMn3kgsI0vYF4rTdjGGGivlXAgWikVDsQAR4CrgLcszwfKa0UEmHaftIFLgQrgFaXUl0qpF5VS\nHYFEYJKlu+pjpdSPAzDGu4AnlFLfA08Cv/NjjLZuBl63/LuX1voIgOVzT79F1cA2Plv/D1jr41gc\nqY9RKXUtcEhrvc2/ITVi+z0MpNeKLdsYA+a1orU+ZInhO4xkXQl8DpyyeQN5EKNFLkQjkrSNd7xX\nAMu01qOBcxjduOFAd+AnwALgr0opFWAx5gB3a60vAe7G0hL3J8s43LXA//k7FnscxaeU+j1QC6zy\nR1xNYqmPUSkVA/weWOTfqBrY+R4G0msFsBtjwLxWLPM6rgMGAX2BjsB0O6cGSq+KCCCStI13tAe1\n1p9aHr+FkSAPAqu1YQtgxqgPHEgx/gJYbTn2f4BfJ9dYTAe+0FoftTw+qpTqA2D57O+u06bxoZT6\nBTADyNBaB8IfStsYL8P4475NKVWG0X3/hVKqd4DEB4H1WrFqGmMgvVauBr7VWldorWsw4hoPdLN0\nl4Pxc/Z3F74IQO0+aWuty4HvlVJDLIfSgN3AGowxJpRSiRiTv/yy4YCTGA8Dky3HrgL2+yG8pubQ\nuOv5bxh/MLF8ftfnETXWKD6l1DQgF7hWa13lt6gaq49Ra71Da91Ta52gtU7ASJBXWH4n/B6fRcC8\nVmw0jTGQXivfAT9RSsVYeiSsr+ePgP+0nBMIrxURgKS4CqCUGgW8iPHH5hvgVowu6JeBUUA18Fut\n9T8DLMbhwDMY3ZMXgKuy2tgAAAI8SURBVHla68/9GGMM8D1wqda60nIsDmNW7ACMP1Y3aa1PBlB8\nXwGRwAnLaZ9oref6Iz5LPM1ibPJ8GTBGa+2XpOjgexhBYL1W7MU4kcB6rTwIzMYYkvkSuA1jDPsN\nINZyLNM6uU8IK0naQgghRJBo993jQgghRLCQpC2EEEIECUnaQgghRJCQpC2EEEIECUnaQgghRJCQ\npC2EEEIECUnaQgghRJAId32KEKIppdTvgDiMoixxwIta6z3+jUoIEeqkpS2EE0qpDkqpR5RSTyul\nllqOjcMok3nG8vlbSdhCCF+QlrYQzmUB0cApjC1SAfYBG4ElQJzW+nn/hCaEaG8kaQthoZQaCTyL\nsUPVUEBh1NSe16QG9Chgm81nIYTwCUnaQgBKqSjgTeAWrfUWpdQfgSigCHhVKfU98E+t9TpgJPAv\nYJLlsxBC+IRsGCIEoJSaAfxca/1zy+ObgGla61/5NzIhhGggE9GEMIwAdtg8vgL4wk+xCCGEXdI9\nLoThBHAVgFIqEbgBGO/XiIQQoglpaQtheB3opJTaCRQAc7TWJ5RS/ZVSHyml7lZKvennGIUQ7ZyM\naQvhhFLqGuBHWutnlFKrtNYZ/o5JCNF+SUtbCOdGAh9a/i3vcIUQfiVJWwjnfgTsU0r1AMr9HYwQ\non2T7nEhhBAiSEhLWwghhAgSkrSFEEKIICFJWwghhAgSkrSFEEKIICFJWwghhAgSkrSFEEKIICFJ\nWwghhAgSkrSFEEKIICFJWwghhAgS/x/p4OhkhHAfCAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10fccee48>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "def create_plot(S, angles=(0,0,0)):\n",
    "\n",
    "    strikes = np.random.rand(10000) * 360\n",
    "    dips =  np.random.rand(10000) * 90\n",
    "\n",
    "    out = np.array([ compute_stress_components_on_fault(S, strikes[i], dips[i], angles) for i in range(len(strikes)) ])\n",
    "\n",
    "    normal = out[:,0]\n",
    "    shear_slip = out[:,1]\n",
    "    shear_dip = out[:,2]\n",
    "    shear_mag = out[:,3]\n",
    "\n",
    "\n",
    "    vert_ss = np.logical_and(dips > 85, np.abs(shear_slip / shear_dip) > 5)\n",
    "    oblique_slip = np.logical_and(0.75 < np.abs(shear_slip / shear_dip), np.abs(shear_slip / shear_dip) < 1.5)\n",
    "    other = np.logical_not(np.logical_or(vert_ss, oblique_slip))\n",
    "\n",
    "    plt.figure()\n",
    "    plt.gca().set_aspect('equal')\n",
    "    plt.plot(normal[other], shear_mag[other], 'k.',\n",
    "             normal[oblique_slip], shear_mag[oblique_slip], 'go',\n",
    "             normal[vert_ss], shear_mag[vert_ss], 'ro'),\n",
    "    plt.xlabel('$\\sigma_n^{eff}$')\n",
    "    plt.ylabel('${\\tau}$');\n",
    "    plt.legend(['Other','Oblique-slip', 'Near vertical'],bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.);\n",
    "    \n",
    "    \n",
    "create_plot(S, angles=(-70,0,0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plot above sort-of speaks for itself, there is no reason to describe in words the answers to the questions with a properly made figure."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem 8 (10 points)\n",
    "\n",
    "Repeat Problem 7 but make $S_{hmin}=S_{Hmax}=45$ MPa (i.e. stress isotropic in the horizontal plane. How has the Mohr diagram changed? Is this still a 3D problem?\n",
    "\n",
    "**Solution**\n",
    "\n",
    "We can use our code above to quickly solve this problem by simply plugging in the new numbers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfMAAADfCAYAAADiFjqyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XtcVGX+B/DPd7gjaiB4w3uJyEWz\nsNZNs5eY2qppua4lqNWWKdav0my1i9rNdUu7uIpm5qXAslVXS61MylprK7ESFdBcxRQv4QUVQW7z\n/P6YgQDnxmXOmWE+79drXsOcmTPny3icL89znuf7iFIKRERE5L4MegdARERE9cNkTkRE5OaYzImI\niNwckzkREZGbYzInIiJyc0zmREREbo7JnIiIyM0xmRMREbk5JnMiIiI3x2RORETk5rz1DsARoaGh\nqlOnTnqHQUTkVnbv3n1GKRVWj/1bent7LwcQAzb+9GYEsK+srOzBG2+88beaT7pFMu/UqRPS09P1\nDoOIyK2IyNH67O/t7b28devW3cPCws4bDAYu5KEjo9EoeXl5UadOnVoO4M6az/MvLSIisiYmLCzs\nIhO5/gwGgwoLC7sAUy/J1c8768AiskJEfhORfVW2hYjI5yLyi/k+2FnHJyKiejMwkbsO87+Fxbzt\nzJb5KgBDamybASBNKdUVQJr5MZHrGDgQSqTytidMkHONVNtWm5tRBBd9Tfc51wjmTorW+zckcjv/\n+9//fOLj46/t2LFjTPv27WPuv//+9leuXJFvv/02YO3atc0rXjd16tS2s2bNaqVnrHpxWjJXSn0N\n4FyNzSMArDb/vBrASGcdn6iqpC1J2BPmQAJOS4MAlbceZ4BOF1BtW21uBgDNSk33nS4AM9/KtBtD\nuQjW3MpOKyIAMBqNGDly5HV33nln/tGjR/cdOXJk3+XLlw2PPfZYeHp6euCWLVua238Xx5SVlTXU\nW2lO62vmrZRSJwHAfN9S4+NTI5e0JQnyvGBRb0GZ4fcEuXjYEvQ441gCrqrm4/py5PheAO79Tz6K\nvKsn+bxAwdhRgoHvDmzgqIgazvbt25vMnDmz9fbt25s0xPt9/PHHTf38/IyPPfbYWQDw9vbG0qVL\nj73//vuhc+bMaffxxx8HR0ZGRr399tvBAJCVlRVw0003dWvXrl3sSy+9VJljkpOTQ2JjY7tHRkZG\njR07tmNF4g4MDOz1+OOPt+3Ro0dkWlpaUEPErAeXHQAnIhNFJF1E0vPy8vQOh1zQ9PvCcdG3esJb\nPGwJjHOAKemAt7KdqF2ZAAgorx57WBGQugH4fEJatd+52FuA1FR9AyaCKZEPGzYs4pVXXgkfNmxY\nREMk9L179wb07NmzsOq2kJAQY3h4eMnUqVNPDh8+/Hx2dnbmQw89dB4ADh065P/VV18d3LVrV9b8\n+fPbFhcXy48//ui/bt26kPT09Ozs7OxMg8Ggli5d2gIAioqKDDExMUUZGRnZgwcPLqhvvHrROpmf\nFpE2AGC+v2quXAWl1DKlVJxSKi4srM7TJKmRSN2bis+7VE/cr6w+gWaljrWwGwtLv6dfOaASE6t9\nNojmtXnSXlpaWtPS0lKD0WhEWVmZIS0trWl931MpBRG5ahCeeftVrx80aFB+QECAatOmTVlISEjp\n8ePHvT/99NOm+/btC+zZs2f3yMjIqJ07dzY7fPiwHwB4eXnhvvvuO1/fOPWmdTL/CMAE888TAGzS\n+PjkRqaMDa4cfDa2RyIGHtE+YbvLMN6an4vKzKzWPb/iSXbNk/PFx8df8vHxMXp5ecHb29sYHx9/\nqb7vGRsbW/Tzzz9Xa+GfO3fOcOrUKV8vL6+r/ov6+flVbvPy8kJZWZkopWT06NFns7OzM7OzszNz\ncnL2vfbaaycAwNfX1+jt7RYlV2xy5tS09wH8F0A3ETkuIn8FMA/A7SLyC4DbzY+JAABzJ0UjL/D3\n1uWi9/OvGnymJQXgjB+Q09z0c11uRujzB0HN7vn7F6RVjq7fMZStdnKOgQMHXt68efPB6dOn527e\nvPngwIEDL9f3Pe+8885LV65cMSxatKgFYBqklpSU1H706NFnWrduXVpQUGA3jw0ZMuTi5s2bg3Nz\nc70B4PTp014HDx70rW9srsSZo9nvVUq1UUr5KKXaKaXeUUqdVUrFK6W6mu9rjnYnD7MrOrgyec98\nKxNhRdokb3tJuEyAxXHA7PWT0SlfQVTdbgalsHJaPC761O4PgPIG/uWrjq7vvzUTn3URyPOCoLlB\nSN3L6+3UcAYOHHj573//+6mGSOQAYDAYsHHjxkMbNmwI7tixY0znzp1j/Pz8jAsXLsy94447Lh08\neDCg6gA4S2688cYrzz77bG58fHxERERE1IABAyKOHTvm0xDxuQpRyvU7EuPi4hTLuTYe0+8Lx/Pv\nnUCA0fTYWUnb2pl9yReYOeoaLF7jOpfJkrYkYUn6ksrH92YAb34ChBZVf11DfVY1PxujAEtuBIJX\npCAhNqGBjkJ6E5HdSqm4uu6/Z8+enJ49e55pyJiofvbs2RPas2fPTjW3u+xodmpc1twaDGOVgWuB\nxoZtfdds3V70ASaM8sKajJSrWsvNipVLJXIASB6aDDVbVd7WrFcIKzTFe/vqeBjmABmh1X/H+rhq\nOpwyzQC4t0ciFvUW+L/kzxY7kRthMienSd2bivF/9kahl+De/+TDgIZJ4Ja6prd3RrXE3axE4d11\nZY2ilbl9/Hao2Qo986r8UZKSggtB3ld9DvVR0Q0/JR0oeq4YY3skosCX096I3AGTOTW4t//gi3KD\naQT66vXlla3wurJ0LbvJiwHVWtu3H1aNInE7LCEBzS+VVv7+azJSsOwmL5RJ/ZN71RZ7UOnv0962\ndREkbUlqqN+AiBoQkzk1mOjF0TjaVPDg96XwUvVrhVcko3KYkvftq+MhSsHbqPDILoXCZwvtvINn\nSYhNwMPfl8HbaEruUzZPRuLdwEXf+rfaK/4dbz8CU1EeEWyMD2+YwImoQTCZU70MfHcg5HnBZ10E\n+x7JRIeC+idwBWBbZ2DK5snwUqbkvX389oYL2gMkD01G6nrT+ABRCt/8fTKONjeNlDfW8T2rjogf\n8cUJKBHsimYNeSJXwGROdRI8LxjyvCDtSBo+XQ0MOlL7lnjN6VibBrSt7DYefFgheWiyM0L3SH1n\nJKNjvmmq3Kpp8cgLaJiu+LjMfCgRFLGkLJGumMypVnZ09YUSwbmZ+TDOQbVE7qiqrW+vOYI1GSkw\nKIWRabnOCZqqeWD+9sqR8qv6BFQWtqlLYq9I6gHmkrIFzQKY1KnBWVsCdeHChS3Gjx/fwdI+gYGB\nvQAgJyfHZ8iQIV20jdik6pKsjz/+eNuNGzfWu7ytNUzm5JBxo7xQLIL+h0qrDZAadMSx/SuSRV4A\nsHJafGXr2zjb6FkD11zM/d8WwmDuDXlqQttqxW1qSwAEXboClZiIQi+u3e6Jlu5aGtJ2QdtYw/OG\nG9suaBu7dNfSkPq+p60lUB3Zv1OnTqWffvrp4frGUV9vvPHGiZEjR9a7vK01TOZkU9KWJIwdJVi+\nyQg/1H6JUAWgFEDC3aZr4GGFCg/M5/VvV/Tqqlw0KzEl9uU3+9S5xS4AAo3AU29lYvp9HCjnKZbu\nWhryxLYnOp4sOOmroHCy4KTvE9ue6FjfhG5tCdS1a9eGFhYWGnJzc3369evXtVOnTjHTpk1rU3P/\nAwcO+Hbt2jUaAAoKCmTYsGFdIiIiooYOHdqlR48ekV9//XUg8HtLHgBWrlwZPGrUqE4AcOLECe/B\ngwdfGxMT0z0mJqb7tm3bLK4El5SUFH7ttddGR0RERE2cOLFdzedHjRrVaeXKlcEAEB4eHjt58uTw\n2NjY7rGxsd337dvnV5/PCADcv7o8OUX4gnDM/OAEFu5G5ch0WxSqv6YiAWzrDGz652Ss4fVvt/LQ\ndyUATLUCwkYk4vZaXkoBTF8ur6w+AbVacMkHaLYyBUhgL0xj9cLXL4RfKbtSrYF4peyK4YWvXwif\n1HtSnUt3W1sCtU2bNiVlZWWSkZHRZO/evfuDgoKMvXr1ihoxYsSFW2+91eJ0l/nz57cMCAgwHjx4\nMPP7778PuOWWW6LsHf/hhx9uP3Xq1NODBw8u+OWXX3wHDx7c9fDhw/urvub06dNeW7duDT58+PA+\ng8GAM2fOeNl732bNmpXv3bs3a9GiRS0effTR9l9++eUhe/vYwpY5VbNjaDTKDILjT56otia4Pds6\nVx/Q9tV1PhzI1ggkxCZg0GFTkZoCX6l1S73ickyzUsCYmMjV2xqxUwWnLC5cYm27o+wtgdq3b9+L\nrVu3Lg8KClJDhw49v2PHjiBr77Vz586gcePGnQWAm2++uSgiIsLuHNdvvvmm2WOPPdYhMjIyavjw\n4dcVFBR4nT9/vlruDAkJKffz8zPec889HVevXn1NUFCQ3UkjEyZMOAcADz300LmffvrJasyOYjIn\nAKYW2NGmgv5bMysTuCNJvGIg25AJvw9mE6Vw2y8lTo6YNJWQgKBio2m2weTJ1brgHU3uBgD3LUjD\n2FEsPtMYtQ5qbfE/vbXtjrK3BGrNNc0trXHuyPNVtxcVFVU+UEohPT09q2L51N9++y0jODjY2Ldv\n366RkZFRY8aM6ejj44Off/45a9SoUfkbN2685rbbbutq7/cyGH5Pv5b+WKktJnPCot6Ce3sk1nqO\nuIKpXviQCUDK3SkczOYpkpNhqFIzvqJXxhEGAKkbTMVnygyC7DFsqTcWs26dlevv7V+tRerv7W+c\ndeusek1TsbUEamBgoHHnzp3NTp8+7VVQUCBbt269pn///gXW3qtv374FKSkpIQCwa9cu/4MHDwZW\nPNeiRYvSH3/80b+8vBybNm0KrrLPxX/84x8tKx5/++23AQCwc+fOX7KzszPXrl179MKFC4Zz5855\njRkz5sLSpUuPZWVlBcKOd999NwQA3nnnneBevXrVe4U5JnMPlrQlCZ91EUxJR2XddHuqjkpPuBvY\n90UK1GwPK6VKAH6vGf/qnPhaJfSKXh9vBXT7MA17W2q9Uj05w6Tek869Puj1o22C2pQIBG2C2pS8\nPuj1o/W5Xg7YXgIVAOLi4grGjBnTOSYmJnr48OHnrV0vB4Ann3zyt8uXL3tFREREzZ07t3VsbGxl\nEn3++edzR4wYcV2fPn26tWrVqrRi+7Jly479+OOPTSIiIqKuvfba6EWLFoXVfN/8/HyvIUOGdI2I\niIjq169ft5deeumYvd+ruLhYevToEZmcnNxq4cKFdl9vD5dA9UBJW5KQv2JJ5RKbjn6VKphKqz46\nzNQSZwKnqnbOS0K3F5bU6pwCTOdVkQH44Il4znRoYFwC1babbrqp2/z584/Z+gPAGcLDw2PT09Oz\n2rRpU1bbfbkEKgEAHhwTiOdHL0HqBiDMwS/digVOlsQBj+xSbImTRX1nJCOs0DRfvRi1a6kHGoH7\nF6ThrD9b6UR1wWTuIVL3pmJPmODtD4scTuKAqczqQ38JgLdRIWmX6/fikP5eXZULP6WwaUBblKN2\nST2kGCgT4ah30sQPP/xwQOtWOQDk5uburUur3BYmcw+w2DzArceZ2nWplwNYNS0ey9dyhTKqvZFp\nufAyr72eG+zt0Mh3gWl++v0L0vBZF8HAd5nUiRzBZN6IJW1Jwp4wQVIdBrhltTTAS7FaGzWAhASE\nnyvFlM2TkRHqWEu9olTw5xPSsKNrvaYpE3kEJvNGyvcFX0TNWeJwa7zqCPU1GSmIOl3u7BDJwyQP\nTUbPPOXwyPeKUe/9D5WixNebC7gQ2cBk3sgkbUnCot6CwjmlmJLuWBnWiiS+7fsUrFnPwW3kXIMP\nKzxy7zUOX08XAL6l5Sgfl8jFW4is0CWZi8gTIrJfRPaJyPsi4q9HHI3N9PvCsXDYEofLsFYUfQmf\n4cMkTppavOY8vJTCsaaOV5HzUsCTb2ciejETuicRkRsfeuihyoVLZs2a1Wrq1Klt9YzJmhkzZrSu\n+rhXr16RdXmfqkunOkrzZC4i4QD+D0CcUioGgBeAe7SOo7GZfl84/rH6BLzhWBKvKMO674sUlMxi\n6VXSR4eLpkpyjna9+xiBfY9k4nwQu91d0tKlIWjbNhYGw41o2zYWS+u/BKqvr6/aunVr8MmTJxt0\nYbDS0lL7L3KQ0WhEeXk5Fi5cWG3Vtp9++im7wQ5ih17d7N4AAkTEG0AggBM6xeH2KqacvbL6hEP/\nmBWtcb/nfTD4MFvjpL/t47dj8GGFg619HBrtLgCCL5dDJSZiYzyXWHUZS5eG4IknOuLkSV8oBZw8\n6YsnnuhY34Tu5eWlxo8fnzd37tyrWqrWlif98ssvA3v16hXZvXv3qF69ekXu2bPHDwAWLlzY4o47\n7ugyYMCA6/r16xdR9b0mT54cPm/evMrqblOnTm07e/bsVgDw3HPPtYqJiekeERER9cQTT7QFTEur\ndunSJToxMbFDdHR01JgxYzoVFxcbIiMjo+68887OQPVlVZ999tlWERERUd26dYtKSkoKB4AFCxaE\nxsTEdO/WrVvU4MGDr7106VKdc7LmyVwplQtgPoBfAZwEcEEptU3rOBqD8AXhiBng2JSzisIvy2/2\nQc88xdY4uZxuJ0vw1Z+iKhdxsUcAjPjiBDLCWGjGJbzwQjiuVF8CFVeuGPDCC/X+i2v69Om/bdiw\nIeTs2bPVlhatWJ503759Wf/+97//N2nSpE4A0LNnzys//PBDdlZWVubs2bNzn3rqqcpu+h9//DHo\n/fffP/Ldd98drPpeiYmJ59avX1/5h8emTZuCExMTz2/YsKHZoUOH/DMyMrKysrIyf/7558BPPvkk\nCABycnL877///rNZWVmZ69aty/Hz8zNmZ2dnfvTRR0eqvveHH37YbMuWLcG7d+/OPnDgQObs2bNP\nAUBCQsL5ffv2ZR04cCCzW7duRQsXLgyt62ek+XrmIhIMYASAzgDyAfxLRBKVUik1XjcRwEQA6NCh\ng9Zhurzk3oJj6Y6tblZRhvXv97RF7rR6rXlA5FS3bTEtE70rOhhxmfl2z20BEHsG+KyLYPBhFjXS\n1SkrS51a214LISEhxtGjR5+dN29ey4CAgMrFXL755ptmv/zyS0DF44rlSc2LnnTOycnxFxFVWlpa\neSr169fvYqtWra6arnPLLbcUnT171jsnJ8fn5MmT3s2bNy/v2rVryauvvtry66+/bhYVFRUFAIWF\nhYbs7Gz/Ll26lLRp06YkPj7e7iIpn3/+ebPExMQzTZs2NQJAxfF3794dMGvWrPBLly55Xb582at/\n//4X6voZaZ7MAQwEcEQplQcAIrIBwB8BVEvmSqllAJYBptrsWgfpypJ7CyY7MFIdMCXyX4OAjY/G\nI3c854yTe+i9/zwyW3mh+29GhxL6oCPA0pu9MOl7TqnUTevWJTh58urE3bp+S6BWmDlz5ukbbrgh\n6p577qmsFV+xPGlQUFC1HPHggw926N+//6XPP//8fwcOHPAdMGBAt4rnAgMDra41Pnz48PMpKSnB\np06d8hk1atS5imM8/vjjJ6dPn16tRv2BAwd8bb1XVRVrr9c0ceLEzuvWrTvUp0+fooULF7b46quv\nmjryfpbocc38VwB/EJFAMf128QCydIjDLY0dVbtEftYP6HhJYTsTObmZqNPlWDktHnkBjlWOe/gH\nI0oMgp3zuFa6LmbNyoV/9SVQ4e9vxKz6LYFaoVWrVuXDhw8/v2bNmsquaGvLk168eNGrXbt2JQDw\n1ltvOdx1PW7cuHPr168P2bx5c3BiYuJ5ALjjjjsuvvfee6EXLlwwAMCRI0d8cnNzLTaEvb29VXFx\n8VVfz0OGDLn43nvvhVZcEz99+rQXYGrld+jQobS4uFg++OCDeo0t0OOa+fcA1gH4EcBecwzLtI7D\nHS3qLUjZ4Hgiz20KhF5hpwa5rwfmb0dYocL4uw2w1wQSAL4K+MPMJaztrodJk87h9dePok2bEogA\nbdqU4PXXj2JS/ZZAreqZZ545lZ+fX5lIrS1P+re//e3UnDlz2t1www2R5eWO99bExcVduXz5sqFV\nq1YlHTt2LAWAu+++++Lo0aPP9e7dOzIiIiLqrrvuujY/P9/L0v4JCQl53bt3rxwAV+HPf/7zxTvu\nuCP/+uuv7x4ZGRn14osvtgaAGTNmnLjpppu69+vXL6Jr165X6vCRVOISqG5gxZMDMTw5zaGlJSum\nnX39p6jK649EjcF9o32xcl2pQ3/MGgHMezgKTy/17P8DXAK18eESqG7qwTGBGPt6mkMrnVUMdDMo\nxUROjc6qf5lGuzvS/DAAmPlWJr6JDHR2WEQugcnchfm+4Iu/f1wEfzv9iwqmlkiyeb1xosbqti37\nISkpuOTr2HX0Px4owt6WnLpGjR+TuYt6cEwgcueVIrTI9uuMMNVVH7Q6HlOYyMkTJCSgabHjRWZi\n8oBjzZjQqXFjMndB40Z5YdG6Irtd6xWt8aFzUjhanTxOt5Ml+LZbgEMJvd0lYG9LQepeloCtJaPR\naORfQi7C/G9hsa+WydzFyPOC1z4x2u1aLxcg8W5TtzpLspKnuiW70KHr6BUt9LARiUzotbMvLy+v\nORO6/oxGo+Tl5TUHsM/S83oUjSEr7hvtiwubgKZ26v8XewH3jwDWrGe3OtFtW/Zj57wk3PjsEviX\nW+/NEgC3HwEu9UoEyvgHsCPKysoePHXq1PJTp07FgI0/vRkB7CsrK3vQ0pOcmuYiEkYJVmwC/OxM\niSwTYPxdTOREluxtKYjJs315SsHUs+X9XgqQ0LiTen2nppH7YMvcBXzWRZByxP7UM7bIiWyL/U3h\nWDNBu0u2W+jeCii9bzx8gEaf0MkzsNtEZ3tbCgY5kMjLBZgyKoCJnMiO9hcVcttfY/c6uk+ZEQUP\njNMkJiJnYzLX0fT7wu12CQKmCyX33WXA8rWFWoRF5Pba/XreoZHuTUoUl1ClRoHd7DoJfCkQZ1OK\nHErkyXHAe+u5GhRRbdySXYiL3oKmdgbFxZ4BjjYVdLzEXi9yX2yZ62DcKC9ceK4I/jbyswKQFwBM\nuNvAqm5EddSsTOF4U9vV4gRAhwJgfyu20Ml9MZlrbM2twXh3gxE+sN5aUAC2dQauf7EtW+RE9dT+\nosLKafE2V10TAFG/gbXcyW0xmWtoxZMDcc9/8u1Om8kIBR6Y0ha50xpkGWAij/fA/O1Y3cf2NfSK\nWu5v/8FXq7CIGgyTuYbueT3NoQ+896M+TOREDez+bwuxN9R+l/uD35diythgrcIiahBM5hpZ1FsQ\n4MDqZ19ea0DJrBJNYiLyND3yFH4Nsp/QX12bz7Kv5FaYzDUwdpRgSrr9qlTHgoABh3iNnMiZOl5S\nyGxpO6EHGIGLf+UcdHIfTOZONnaUYOUm+4l8cRzQgVNjiDQRfVrZnIcuACbtUsi8PlzLsIjqjMnc\niYLnBePNT2zXW1cwrUfO6WdE2rol23YRJgHQfc8JHGjDAXHk+pjMnSR6cTTu2JWP0CLrr6mYgsYS\nrUT6OHGz7eVTBUDEqVLs6MqETq6NydwJUvemoucXmXj7Y9tzyRfHAYMPM5ET6SX8u/1Ij7Jdx10A\n9D9UioHvDtQqLKJa0yWZi8g1IrJORLJFJEtE+ugRh7NsmZOI9zYATWysS34mAAhekaJdUERkUe/9\n57FpQFu7ddxbbkrjCHdyWXq1zN8E8KlSKhJATwBZOsXR4MaOEry7AfCy8ZorBuDt8VFIiOXSi0Su\nYGRaLrJ6Wk/oAiBlA1A26WEtwyJymObJXESaAbgVwDsAoJQqUUrlax2HM/i+4It3NtlevaZMgGlj\nrsHTS/drFhcR2Rf1c67NojIGAOO/vYxd0SwoQ65Hj5Z5FwB5AFaKyE8islxEmugQR4NK3ZuKj1eV\n2l08ZfxdwOI15zWLi4gct/eLFGzrbD2hC4C4zHwOiCOXo0cy9wZwA4AlSqleAC4DmFHzRSIyUUTS\nRSQ9Ly9P6xhr7fwDiRh0xPZ88jMBHLlO5MoSYhMwdVoUjja3/pqKAXErnuSAOHIdeiTz4wCOK6W+\nNz9eB1Nyr0YptUwpFaeUigsLC9M0wNryet4Lkx2o8Pb2+CitQiKiOto/ZT/W3BNld5W14clpWoVE\nZJfmyVwpdQrAMRHpZt4UDyBT6zgaSureVIzJMNr8ICvmk/M6OZF7eHrpfhgmT7Y5wj20CNgYzwpx\n5Br0Gs3+KIBUEckAcD2AuTrFUW9b5iTi3X/bnk/+eWfOJydyO8nJdq+fj/jiBJCUpGVURBbpksyV\nUj+bu9B7KKVGKqXcckTYlLHBWLkJ8Lbyv72iRT6IiZzILW3652QU2fiWFADGJUs0i4fIGlaAq6OB\n7w7E/A/zbdZdv+QLnNnEwjBE7ip5aDJmj7NdUEYALshCumMyr6Ppc9JsTkO77AP8350+LAxD5OZe\nXZWLc4HWh7dWLMiycx6720k/TOZ1MG6Ul81paArAQ8OBVf8q0TIsInKSFsveQ5mN6SoCoOdsdreT\nfpjMayl6cTRWbTByPjmRJ0lIgPd7KTa724NKwNY56YbJvJZeW5Bpdxrax0nxWoVDRFpJSLBZTEIA\ndH9xqWbhEFVlN5mLyC4ReUdEHheRASLi2hVcnGjupGi73evbOgMPzN+uZVhEpBGZZHvueUihwo6h\n0ZrFQ1TBkZb5CAD/AuALYBKAHBE56tSoXNS0tzJtdq8DwANT2moSCxHpIDkZEh9vc+75rVszWeqV\nNGc3mSulTiilPlVKvaKU+guAOADLnR+aa0nuLbC1tIICkBEK5E7L1SokItLD9u145w++NldXY6lX\n0poj3ewdqj5WSmUB8Kh+pIHvDsQkG7XXFYAzfkDPPA56I/IEActW4EyA9edDi0zfG0RacaSbfa2I\nHBeR/4hIsoi8BiDS2YG5kpH/TLO7iMq2XSwOQ+QpEmITcGCW7evnLTelIXVvqmYxkWcTpRxrTYrI\ndQBiAYQA+EwpddyZgVUVFxen0tPTtTpcNYEvBaLguSKrf/UoAIvjgEd2sVVO5Gku+gmaWSknkdMc\niHmqCQqeLtA2qCpEZLdSKk63AEgzjnSztxCRyQD6AcgFsEbLRK63VzYV2WyVl4CJnMhTrXvU+mC4\njheAO3dfZuucNOFIN/u/AYTPsDEpAAASJUlEQVTBtLLZqwAuiEi2U6NyEeELwpFk51r5lL/YuHBG\nRI3aA/O347yVUq8CIGUDUDTxAW2DIo/kSDJvqpR6AcBppVR/APcCWOncsFzDt3NO2J1TvnxtoZYh\nEZGLCVn2Hi77WH7OAOCv35WwMhw5nSPJ/Ir5vlhEApRS6wH8yYkxuYSMMEGHApsFn7giGhEBCQmY\nONz2uuc3Psu67eRcjiTz+SISAmAtgBUi8iiARr/eX+wZ24lcAVwRjYgAAM0fmIyjza0/71/Ouu3k\nXI4UjVmvlDqnlHoNwFYA7WGqCtdoyfO267wpAIbJk7UJhohcXvLQZHzxoO3KcDc+x9Y5OY/DU9P0\npPXUtD1hgh5WWuYKwJfXGjDgkI3FzInIIxlFbE5j/epPUbhty37N4uHUNM/BVdNq2BgfbjWRA0A5\nwERORBYd/Iv9uu1EzsBkXsOwL22PYH/wz1aGrRKRx4tcux2lXtYv09lbqImorpjMq4heHA0vG1cd\njABW/ctKuSciIgC+q9+zWeY1ewxrtlPDYzKv4rUF1rvAFIAlvPJERPYkJNjsau/2YRqQyqpw1LB0\nS+Yi4iUiP4nIZr1iqGrK2GAMOmJ90FtGKJA5hyPYici+jwa0tZnQi/46QctwyAPo2TJ/DECWjsev\nZv6H+TavZ13/iGn6CRGRPSPTcnHJxvAa/+JyrHiS3e3UcHRJ5iLSDsBQAMv1OH5NcydFw9/GAPVy\nAdoGtdUuICJyexkvWF8iVQDc+0aaluFQI6dXy/wNAE/BNKbMIhGZKCLpIpKel5fn1GDGfpBpcwT7\n0huB3Gm5To2BiBqXvjOSURrgZ/V5/3Jgx9BoDSOixkzzZC4iwwD8ppTabet1SqllSqk4pVRcWFiY\nU2PqcMFGHOASp0RUN75vv2Ozdd73E847p4ahR8v8FgB3ikgOgA8ADBAR3VYsCZ4XbPU5BSCZI9iJ\nqK4SbK/fYGsqLFFtaJ7MlVIzlVLtlFKdANwD4AulVKLWcVT4YJntgW/PjbpGs1iIqPE5F6B3BOQJ\nPHqeefiCcKvT0QDgaHPg/IzzmsZERI1L1izrA+HYMKeGomsyV0rtUEoN0+v4/b89YfU5BeDpeO1i\nIaLGqe+MZJwNtPECFpChBuDRLfNVG23XSl6znn83E1H9ZT832eLUHQOAvIkTtQ6HGiGPTebhC8Lh\nY2VinAKwrbOm4RBRI9Z3RrLVhkOLwkIkbUnSNB5qfDw2mc/8wHoXOwBMnRalUSRE5AmkY0eL239t\nDixJX8KETvXimck8NRVT0m13se+fsl+zcIjIA7z8MhBY/eL5ZZ/fx+Ys271Mh6CosfDIZH7+4Qk2\nK75lhGoZDRF5hIQEYNky5DQ3lb7MCwAKvYGUDcCR14G/7LFRU5rIDo9M5tdctv2fpmceB74RkRMk\nJKDzE0Di3UCzEiCsyPQl3OkCsHITOLKd6szjkvncSbZrIZ9hgQcicrI3PwH8arQp/MoBPPaYLvGQ\n+/O4ZP7ICtuLqrw9ngPfiMh5OjbviNAiK0+ePatpLNR4eFwyb1pq+/mnl3LgGxE5z8vxL1t9jhf4\nqK48Lpnbwi52InK2hFjbi68Q1YVHJfMHx1ivqcgudiLSiq2Gw855nG9OtedRyfzZz4pszi1nFzsR\naeHFu1tY7FIXAJEvLtE6HGoEPCqZd7hg/blyW1meiKgB3fy3N60+16IQrAZHteYxyTx4XjDOWuna\nUgB23sEudiLShr3r5tFz2Dqn2vGYZJ5fnG/1uSIDcNsWdrETkXaKmzexuF0ATE7XNhZyfx6TzAGg\nhZW5nf5WVk8jInIW/8VvWZ2Kxqt+VFsek8w/XW39P8ivzTUNhYjIVKvdhtS9LO1KjvOIZD53UjQG\nHbGczBWAxSPbah0SEZHNFvj3/2BpV3KcRyTzsR9YL+EKAK+uytUsFiKiSpMnW52i9twGlnYlx3lE\nMrc1JY2ISDfJyVafCi1iVzs5ziOSuS2XfPSOgIjIsoc/fljvEMhNaJ7MRaS9iHwpIlkisl9EnHph\naMWTA22ukpbxwmRnHp6IyKZzgZa/oc4EAJdLL2scDbkrPVrmZQCmKaW6A/gDgCki4rSKLQOWp9m8\nXt53hvVuLiIiZ8t6bhKu1PgmrriOfm+G5uGQm9I8mSulTiqlfjT/fAlAFoBwZx2P18uJyJX1nZGM\nB0YCeQG/J3EBEFYEvP0xgFReNyf7dL1mLiKdAPQC8L2zjlFg45o4lzwlIlfwfg/gsu/VU9WalAJ4\n5hk9QiI3o1syF5EgAOsBPK6Uumjh+Ykiki4i6Xl5eXU+TlCp5e0KwPShHP1GRK7BWi+i+vWotoGQ\nW9IlmYuID0yJPFUptcHSa5RSy5RScUqpuLCwsDodZ+e8JJvXy1f9q6RO70tE1JCa+DSxWonyWHMW\ndyX79BjNLgDeAZCllHrNmcdqP2+J1WTOJU+JyFW8NfwtPB0PXK7RWXjZB5gxwFoFd6Lf6dEyvwXA\nOAADRORn8+1PzjhQe2vdVgCW3uiMIxIR1V5CbALe7wE8NBzIaQ4YYbp/aLjpejqRPd5aH1AptRM6\nLwqkADw6DHhEzyCIiGp4vweTN9VNo64AZ+0vBgEQFeq0qe1ERLXWxMfy+ub+4q9xJOSOGnUyt2X/\nlP16h0BEVOmt4W/BAAPuzQCOvA6UzzHdz9vLpjrZ12iTORcoICJ3khCbgK/8Hsbyj4FOF0xfzp0u\nAI98vIeFY8iuRpvMt8xJ1DsEIqJa6bt0KwJr1MbwKi5m4Riyq9Em8zc/0XmUHRFRbf36a+22E5k1\n2mQeWmT9Oc7aJCKX1MTyIDir24nMGm0yt0YBSI7TOwoiIgsuW17y1FhQoHEg5G48LpkDwCO72DYn\nIhekLH83CTiol2zzyGROROSSvLwsbi4X4Jk0DoIj65jMiYhcxcSJV43pqSg//esFDoIj65jMiYhc\nRXIy3v1jE5SJKYmXCbA4zlR+ukPzDnpHRy5M89rsRERknffStxC44X6U4vcJ5z7wwcvxL+sYFbm6\nRtkynzspWu8QiIjqJCE2AU/HPA25YGqeywXB0zFPIyE2Qe/QyIWJsjJ60pXExcWp9PR0h19/0U/Q\nrMTycwqAuMHvTESe7b///S927NiB2267DX369KnTe4jIbqUUJ+N6gEbZzd7USiInInILqano88wz\n6PPrr0CHDsDLLwMJbJmTdY2ym90aBSAjVO8oiIhsSE0FJk4Ejh41zTs/etT0mIutkA0elcwB4K1V\nk/UOgYjIumeeAQoLq28rLORiK2STxyXz5KHJeodARGQdF1uhOvC4ZE5E5NI6WJlPbm07EZjMiYhc\ny8svo9i7+tjkYm9v0yA4IiuYzImIXElCAlL790cOACOAHACp/ftzNDvZxGRORORiur/4IiL9/OAt\ngkg/P3R/8UW9QyIXp8s8cxEZAuBNAF4Aliul5ukRBxGRK+rTpw++/PLLeheNIc+heTIXES8AiwHc\nDuA4gF0i8pFSKlPrWIiIXFWfPn2YxMlhenSz3wTgkFLqsFKqBMAHAEboEAcREVGjoEcyDwdwrMrj\n4+Zt1YjIRBFJF5H0vLw8zYIjIiJyN3okc7Gw7aqVT5RSy5RScUqpuLCwsFod4ExA7bYTERG5Mz2S\n+XEA7as8bgfgREMe4JW/tMWVGr/ZFYNpOxERUWOjRzLfBaCriHQWEV8A9wD4qCEP8OqqXDw3ri1y\nmpvnaTYHnhvXFq+uym3IwxAREbkEzUezK6XKROQRAJ/BNDVthVJqf0Mfp2ri7gTg1YY+ABERkYvQ\nZZ65UmorgK16HJuIiKixYQU4IiIiN8dkTkRE5OZEqatmhbkcEckDcLSOu4cCONOA4TQUxlU7jKt2\nGFftNNa4Oiqlaje3l9ySWyTz+hCRdKVUnN5x1MS4aodx1Q7jqh3GRe6O3exERERujsmciIjIzXlC\nMl+mdwBWMK7aYVy1w7hqh3GRW2v018yJiIgaO09omRMRETVqbp3MRcRLRH4Skc3mx6tE5IiI/Gy+\nXW9lvwki8ov5NkGDuP5TJaYTIrLRyn7lVV7XoPXqze+fIyJ7ze+fbt4WIiKfmz+Lz0Uk2Mq+TvvM\nrMT1qohki0iGiPxbRK5xdF8nxzVHRHKr/Dv9ycq+Q0TkgIgcEpEZGsS1tkpMOSLys6P7NmBc14jI\nOvO/W5aI9HGR88tSXK5wflmKS/fzi9yUUsptbwCmAlgDYLP58SoAf7azTwiAw+b7YPPPwc6Mq8Zz\n6wGMt7JfgZM/rxwAoTW2vQJghvnnGQD+ofVnZiWuQQC8zT//w1Jc1vZ1clxzADxpZz8vAP8D0AWA\nL4A9AKKcGVeN5xcAmKXD57UawIPmn30BXOMi55eluFzh/LIUl+7nF2/ueXPblrmItAMwFMDyWu46\nGMDnSqlzSqnzAD4HMESLuESkKYABACy2zHUyAqYvFZjvR1p4jVM/M0uUUtuUUmXmh9/BtFSuu7gJ\nwCGl1GGlVAmAD2D6nJ1ORATAXwC8r8Xxqhy3GYBbAbwDAEqpEqVUPnQ+v6zFpff5ZePzcoRu5xe5\nLrdN5gDeAPAUTKucVvWyuevsdRHxs7BfOIBjVR4fN29zdlwAcBeANKXURSv7+otIuoh8JyKWvvTq\nSwHYJiK7RWSieVsrpdRJADDft7Swn7M/M0txVfUAgE/quK8z4nrEfI6tsNJtrOfn1Q/AaaXUL3XY\ntz66AMgDsFJMl5iWi0gT6H9+WYurKj3OL1tx6X1+kRtyy2QuIsMA/KaU2l3jqZkAIgH0hqnL7m+W\ndrewrUGG9NuIq8K9sN1i6qBM1Z7GAnhDRK5tiLiquEUpdQOAOwBMEZFbHdzPaZ+ZmdW4ROQZAGUA\nUmu7r5PiWgLgWgDXAzgJU5d2Tbp9XrB/jjnr8/IGcAOAJUqpXgAuw9St7ghnfl4249Lx/LIWlyuc\nX+SG3DKZA7gFwJ0ikgNTF9MAEUlRSp1UJsUAVsLUHVXTcQDtqzxuB+CEM+MCABFpYY5ni7WdlVIn\nzPeHAewA0KuB4qr5/r8B+Lc5ntMi0sYcYxsAv1nY1ZmfmbW4YB4INQxAglLK4peVtX2dFZdS6rRS\nqlwpZQTwtpXj6fV5eQO4G8Da2u7bAI4DOK6U+t78eB1MyUrv88taXHqfXxbjcoXzi9yTWyZzpdRM\npVQ7pVQnAPcA+EIplVjlS0Nguja3z8LunwEYJCLB5i6sQeZtTovL/PRomAbEXbG0rzkeP/PPoTD9\nYZDZEHGZ37OJ+Zo9zN15g2D6fD4CUDF6eAKATRZ2d9pnZi0uERkCU8/KnUqpwlr+Ts6Mq02Vl91l\n5Xi7AHQVkc4i4gvTudAgsxPs/M4DAWQrpY7XYd96UUqdAnBMRLqZN8XDdP7qen5Zi0vv88tGXLqe\nX+TG9B6BV98bgNvw+2j2LwDshek/QAqAIPP2OADLq+zzAIBD5tv9zo7L/HgHgCE1XlMZF4A/mmPf\nY77/awPH08X83nsA7AfwjHl7CwBpAH4x34do+ZnZiOsQTNcFfzbflpq3twWw1da+To7rPfO/TwZM\nX6BtasZlfvwnAAdhGnXs9LjMz60CMKnG6zX5vMzvfz2AdPNnsxGmkem6nl824tL1/LIRl67nF2/u\ne2MFOCIiIjfnlt3sRERE9DsmcyIiIjfHZE5EROTmmMyJiIjcHJM5ERGRm2MyJyIicnNM5kRERG7O\nW+8AiBoTEZkJU6GUs+b75UqpbH2jIqLGji1zojoQER8RmSsib4hIsnnbzTAtdHLJfH+EiZyItMCW\nOVHdTAQQACAfprKfgKm85g4ACwG0UEot1ic0IvI0TOZEdohITwD/BBAK0xK7AmAFgCRlWqGvwvUw\n1fGuuCci0gSTOZENIuIP05Ki45VSP4jIiwD8AXwNYJWIHINpdbxPAfQE8B8A/cz3RESa4EIrRDaI\nyDAAY5VSY82PR8O0+t1f9Y2MiOh3HABHZFsMTEtSVrgBwI86xUJEZBG72YlsOwtgAACISASAu2Fa\ne56IyGWwZU5k2/sAgkRkH4BlAO5VSp0VkXYi8qWIPCEia3WOkYg8HK+ZE9WBiAwFcJ1S6k0RSVVK\nJegdExF5LrbMieqmJ4DPzD/zL2Ii0hWTOVHdXAfgoIiEAjildzBE5NnYzU5EROTm2DInIiJyc0zm\nREREbo7JnIiIyM0xmRMREbk5JnMiIiI3x2RORETk5pjMiYiI3ByTORERkZtjMiciInJz/w+Tq96/\nXZXVpAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10fc96470>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "S = np.diag([65, 45, 45])\n",
    "create_plot(S, angles=(-70,0,0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This problem is clearly 2D as it can be fully characterized by the two principle stresses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(39.999779996773242,\n",
       " 8.6601270189221928,\n",
       " 6.1230992826130353e-16,\n",
       " 8.6601270189221928)"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compute_stress_components_on_fault(np.array([[30, -8.66, 0],[-8.66, 40, 0],[0,0,30]]), 60, 90, (0,0,0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3875.0000000000005,\n",
       " -433.01270189221918,\n",
       " -649.51905283832912,\n",
       " 780.62474979979982)"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compute_stress_components_on_fault(np.diag([4000, 3000, 5000]), 45, 60, (0,0,0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(4125.0, -433.01270189221964, 649.51905283832934, 780.62474979980027)"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compute_stress_components_on_fault(np.diag([5000, 4000, 3000]), 225, 60, (0,0,0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1441.2168789736368, 488.4593858091485, 160.58981080223731, 514.18057034328797)"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compute_stress_components_on_fault(np.array([[2100, -520, 0],[-520, 1500, 0],[0,0,1000]]), 120, 70, (0,0,0))"
   ]
  }
 ],
 "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}