{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem 1 (60 points)\n",
    "\n",
    "Given the reservoir conditions below, and a desire to keep beakout widths as small as possible for wellbore stability, which direction and type of well (e.g. vertical, horizotonal, $\\phi = 20^\\circ$, etc.) would be the safest to drill.  Support your recommendation with a figure.  Preferably one that quantatively shows the expected wellbore breakout widths for a given drilling trajetory.\n",
    "\n",
    "You can assume a Mohr-Coloumb failure criterion with an unconfined compressive strength of the rock of $145$ MPa and an internal friction of $\\mu_I = 1$.  The Poisson ratio is $\\nu = 0.2$.\n",
    "\n",
    "$S_{Hmax} = 115$ MPa in the direction N$30^\\circ$E\n",
    "\n",
    "$S_{hmin} = 90$ MPa\n",
    "\n",
    "$S_{v} = 70$ MPa\n",
    "\n",
    "$P_p = P_m = 33$ MPa"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution**\n",
    "\n",
    "First we will define the functions we need for all of the operations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "\n",
    "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))\n",
    "\n",
    "\n",
    "def compute_Sb(S, angles_G=(0,0,0), angles_B=(0,0)):\n",
    "    \n",
    "    delta, phi = np.radians(angles_B)\n",
    "    \n",
    "    Rb = np.array([[ np.cos(delta) * np.cos(phi),  np.sin(delta) * np.cos(phi), -np.sin(phi)],\n",
    "                   [ -np.sin(delta), np.cos(delta),  0.0],\n",
    "                   [ np.cos(delta) * np.sin(phi), np.sin(delta) * np.sin(phi), np.cos(phi)]])\n",
    "    \n",
    "    Sg = compute_Sg(S, angles_G)\n",
    "    \n",
    "    return np.dot(Rb, np.dot(Sg,Rb.T))\n",
    "\n",
    "\n",
    "def compute_wellbore_stress(S, nu, theta, DP):\n",
    "    \n",
    "    theta = np.radians(theta)\n",
    "    \n",
    "    sZZ = (S[2,2] -  2. * nu * (S[0,0] - S[1,1]) * np.cos(2 * theta) - \n",
    "           4. * nu * S[0,1] * np.sin(2 * theta))\n",
    "    \n",
    "    stt = (S[0,0] + S[1,1] - 2. * (S[0,0] - S[1,1]) * np.cos(2 * theta) - \n",
    "           4 * S[0,1] * np.sin(2. * theta) - DP)\n",
    "    \n",
    "    ttz = 2. * (S[1,2] * np.cos(theta) - S[0,2] * np.sin(theta))\n",
    "    \n",
    "    srr = DP\n",
    "    \n",
    "    return (sZZ, stt, ttz, srr)\n",
    "\n",
    "\n",
    "def compute_max_tangent_stress(S, nu, theta, DP):\n",
    "    \n",
    "    sZZ, stt, ttz, srr = compute_wellbore_stress(S, nu, theta, DP)\n",
    "    \n",
    "    return 0.5 * (sZZ + stt + np.sqrt((sZZ - stt) ** 2. + 4.0 * ttz ** 2.))\n",
    "\n",
    "\n",
    "def compute_min_tangent_stress(S, nu, theta, DP):\n",
    "    \n",
    "    sZZ, stt, ttz, srr = compute_wellbore_stress(S, nu, theta, DP)\n",
    "    \n",
    "    return 0.5 * (sZZ + stt - np.sqrt((sZZ - stt) ** 2. + 4.0 * ttz ** 2.))\n",
    "\n",
    "\n",
    "def compute_breakout_width(S, Pp, Pm, nu, C0, mu, angles_G=(0,0,0), angles_B=(0,0)):\n",
    "    \n",
    "    Sb = compute_Sb(S, angles_G, angles_B)\n",
    "    \n",
    "    Sb_eff = Sb - Pp * np.eye(3)\n",
    "    \n",
    "    theta = np.linspace(0, 360, num=90)\n",
    "    \n",
    "    smax = np.array([ compute_max_tangent_stress(Sb_eff, nu, i, (Pm-Pp)) for i in theta])\n",
    "    smin = np.array([ compute_min_tangent_stress(Sb_eff, nu, i, (Pm-Pp)) for i in theta])\n",
    "    DP_arr = np.ones_like(smin) * (Pm-Pp)\n",
    "    \n",
    "    sigma1 = (np.array([smax, smin, DP_arr]).T).max(axis=1)\n",
    "    sigma3 = (np.array([smax, smin, DP_arr]).T).min(axis=1)\n",
    "    \n",
    "    breakout_bool_array = C0 < (sigma1) - ( np.sqrt(mu ** 2 + 1) + mu ) ** 2. * (sigma3)\n",
    "    \n",
    "    return np.round(breakout_bool_array.sum() * 2.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "S = np.array([[115, 0, 0],[0,90,0],[0,0,70]])\n",
    "\n",
    "delta = np.linspace(0, 360, num=50)\n",
    "phi = np.linspace(0, 90, num=50)\n",
    "\n",
    "res = np.array([ (np.sin(np.radians(d))*np.sin(np.radians(p)), \n",
    "                  np.cos(np.radians(d))*np.sin(np.radians(p)), \n",
    "                  compute_breakout_width(S, 33, 33, 0.2, 145, 1.0, angles_G=(30,0,0), angles_B=(d,p))) \n",
    "                  for d in delta for p in phi])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATwAAAEICAYAAADC7ki9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJztnXu0HWV99z9fQyFqURIS0nCJCSUV\notaop0GrL1ouCrYvwaotl9rYQvOyKK3F1xZoq1gUQetaoa3XqEi0vICFWqKCNBIuiyrISRshFzDh\nVkJiCAGsigkSfu8f82wymTOz9+w9s/eemf37rLXX2fPMMzPP7LPP5/ye68jMcBzHGQVeMOwCOI7j\nDAoXnuM4I4MLz3GckcGF5zjOyODCcxxnZHDhOY4zMrjwnI5IukXSGcMuRxaSTpP07232ty2/pMsl\nfbQ/pXOqhAuvhkh6SNLPJf1U0pOSviXpkGGXq1vCfRxb9DxmdoWZvTXnNd8r6fai13TqiQuvvvxv\nM/tlYCawFfinrIySJg2sVI5TYVx4NcfMdgDXAPNaaaGK9llJ10v6GfBbkvaR9ElJ/y1pq6TPSXph\nyD9F0jclbQsR4zclHZx2PUkzJd0t6QNh+0BJyyU9IWmjpD9JlOOjse23SNoU3n8VmAV8I0Sqf5Vy\nrVslvTO8f5Mkk/T2sH2spNXh/R5Rm6TjJN0r6ceSPgUopB8BfA54Q7jmU7HLTQmR8k8k3SnpV7v4\nNTg1wYVXcyS9CPh94I7ErlOBi4B9gduBjwO/BswHDgMOAj4U8r4A+DLwMiIJ/Rz4VMq1ZgO3Ap8y\ns0+G5CuBTcCBwLuAj0k6plO5zew9wH8TIlUz+0RKtluBt4T3RwEPAG+Obd+aUsZpwLXA3wLTgPuB\nN4ZrrgfOBL4Xrrlf7NBTgL8DpgAbiT47p2G48OrLv4UI5X+A44C/T+y/zsz+w8yeA3YCfwKcY2ZP\nmNlPgI8BJwOY2XYzu9bMng77LmK3WFrMA24BLjCzpQCh3fBNwLlmtsPMVgNfBN5T0j3eyp6Cuzi2\n/WZShAe8HVhnZteY2S+AS4Ef5bjWv5rZ983sWeAKon8MTsNw4dWXk0KEsg9wNnCrpF+J7X8k9n46\n8CJglaSngii/HdKR9CJJn5f0sKT/AW4D9ku0/Z0GPEpUfW5xINASaIuHiaLHMvge8GuSZhAJ6CvA\nISGKWxDKmeRAYvdu0eoYj6TkSxKX4tPAL/daaKe6uPBqjpntMrN/BXYRRVvP74q9f5yomvoKM9sv\nvF4aOj0A/i/wcuBIM3sJUTQFoe0r8OFwnv8XE+FmYKqkfWP5ZhGJEeBnRKJtERdysoxp9/Y0sAp4\nH7DGzJ4Bvgu8H7jfzB5POWwL8HyPtSTFtztd02k2Lryao4iFRG1P69PyhGrtF4Alkg4Ixx0k6W0h\ny75EQnxK0lTggpTT/AJ4N/Bi4KuSXmBmjxAJ6GJJkyX9OnA6UZUQYDXwdklTQ/T5F4lzbgUO7XCL\ntxIi2LB9S2I7ybeAV0j6XUl7AX/OnqLdChwsae8O13UaiAuvvnxD0k+J2vAuAhaZ2do2+c8laoy/\nI1Rbv0MU1UHUzvVCogjuDqLq7gRChPW7wAHAZZJeQNTYP5so2vs6URvfinDIV4EfAA8B/w5cnTjl\nxcDfhmr2BzLKfSuRkG/L2E6W8XEiMV8CbAfmAv8Ry7ISWAv8SFJahOg0GPkCoI7jjAoe4TmOMzKU\nIjxJl0l6TNKajP2S9I9hYOrdkl4b27dI0obwWlRGeRzHaR6S3idpjaS1kv4ipE2VtCL4Y4WkKe3O\nUVaEdzlwfJv9JxC1pcwFFgOfbRWWqIH8SKJhBhd0KrDjOKOHpFcSjSVdALwa+B1Jc4HzgJvMbC5w\nU9jOpBThmdltwBNtsiwEvmIRdxCN8ZoJvA1YEQbDPgmsoL04HccZTY4A7giD458l6rx6B5FbloU8\ny4CT2p1kr74WcTcHsefgz00hLSt9ApIWE0WHvPjFL37d4Ycf3p+SOj1xzxPtJzO8YIcy973iwBll\nF8cpyKpVqx43s+lFzvGGN7/QnnryuVx5773nmbXAjljS0taMnsAa4CJJ+xMNoXo7MA7MMLMtAGa2\npTXsKotBCS/t225t0icmRje/FGBsbMzGx8fLK53TljlXXNwxz4E5zjP53smp6T/LyL//ul3Pv7/9\n2qxRK04/kPRw0XM89eRzXP6NPN8MeP3sh3aY2VjWfjNbL+njRLXAnxINd3q22zINSnib2HO0+8FE\n47Y2sXtyeCv9lgGVyUmQR2zD4k3v/OSEtH1XPcoNDy0ZQmmcYWBmXwK+BCDpY0T+2CppZojuZgKP\ntTvHoIS3HDhb0lVEHRQ/DgW8kWh1jVZHxVuB8wdUppFm0fdPT0k9bODlKMoJs8+ZkOYSbCaSDjCz\nxyTNIhoA/wZgDrCIaKD5IuC6ducoRXiSriSK1KaF9c4uAH4JwMw+B1xPVOfeSDQx+4/CvickfQS4\nK5zqQjNr1/nhFCBdcrs5au7G59/ftmG48otXZ9PYd9WjmfviEnT5NYprQxveL4A/NbMnJV0CfE3S\n6UTLjb273QlKEZ6ZndJhvwF/mrHvMuCyMsrhTKST5LI4au7G0qW34/Adme14SbbPm9RRenlw+TUH\nM/tfKWnbgY7rL7YYVJXWGSC9Sq7puPwcF15D6IfkWlXcYVdv+4HLbzRx4dWcQURz/ajedqKM6mxe\nWvJz8TUfF14NGUaVtawOjW7a8bJo12FRBI/6mo8Lr0aMUttcWZ0WveJRXzNx4dWApomujChvULj4\nmoULr8JUVXTDaNMbNi6+ZuDCqyBVFV2cotKrU5QXx8VXb1x4FaIOoovT5GErnXDx1RNf4r0i1E12\nZbDj8B2dM1WctLm8TnXxCG/INEF0/WrTS+up7deQlCJ4tFcfPMIbEou+f3ojZNciPk6vG7qN8n7y\nutT1YSvBCbPP8Yiv4rjwhkCTRBfnqLkbexZfk3DpVRev0g6QpoquKHXtsW2HV3OriUd4A2KUZOdR\n3m482qsWLrwBMEqya9Ft9TarLW/7vEl7bFex06IT3rZXHVx4faRpHRP9ZsfhOzp2YlS506ITLr3h\n48LrEy66CO/I2BOX3nApRXiSjpd0n6SNkiY8+VvSEkmrw+uHkp6K7dsV27e8jPIMG5ed0w6X3vAo\nLDxJk4BPAycA84BTJM2L5zGzc8xsvpnNB/4J+NfY7p+39pnZiUXLM2xcdsVJVmu3z5u0R1tenau1\nLVx6w6GMCG8BsNHMHjCzZ4CrgIVt8p8CXFnCdSuFt9e1p9tqbVpbXlJ6dRefd2YMnjKEdxDwSGx7\nU0ibgKSXET1HcmUsebKkcUl3SDqphPIMHBddPsqWHni053RHGcJTSppl5D0ZuMbM4hMkZ5nZGHAq\ncKmkX029iLQ4iHF827ZtxUpcImXL7qwZKzlrxsrOGWtKGR0YLj2nV8oQ3ibgkNj2wcDmjLwnk6jO\nmtnm8PMB4BbgNWkHmtlSMxszs7Hp06cXLXMplCm7posuTj/G57n0nDyUIby7gLmS5kjam0hqE3pb\nJb0cmAJ8L5Y2RdI+4f004I3AuhLK1HfKkF1LckVFV0dZdjNcJWt8nkvP6ZbCc2nN7FlJZwM3ApOA\ny8xsraQLgXEza8nvFOAqM4tXd48APi/pOSL5XmJmlRdeXtmdNWMln9l6dOa+LLKO6XR8Ky1+fNZ1\n8lxjEHSziGjanNvkElI/ed1BtZyNEeeE2ef4HNw+UcriAWZ2PXB9Iu1Die0Ppxz3XeBVZZRhUOSR\nXVwyadIrIru8kVxRoQ6avOLLI70m4NLrD75aShd0K7tu97dElMzzma1Hd1VlrZvs4uRZTNSlN3qE\nJrGrY0mHAh8C9gP+BGj1ZP51CMBSceHlpJPs8lQdO4koa38Z7XNVF12cXqXXognV2hYuvQgzuw+Y\nD89PdngU+DrwR8ASM/tknvO48HJQVHZ5RNfvTodO56+aEHt5QFA8ymt1YDRBfC69CRwD3G9mD0tp\no+Ky8cUDOjAI2Tm9kdZ728SeW/De2wTJ4W1nS7pb0mWSprQ7UHt2mtaDsbExGx8f7/t1Rkl27co8\n7Ogvb5SXrOLG2/TqHum1xH37tR8YyPUkrQoTAnpm2hHT7LeXtZtlupuvHHnZw8DjsaSlZrY0pVx7\nE43zfYWZbZU0IxxnwEeAmWb2x1nX8Qgvg37KLs/+QVPlsXx5x+y1i/bqHOnFy/6md+Zqqqojj7cm\nFoTXBNkFTgD+08y2ApjZVjPbZWbPAV8gmtufibfhdUHedrCqiqPu5GnXa0mvFe21pLf/ul217Myo\ns6j7xB6Lj0iaaWZbwuY7gDXtDvYIL4W06G6UZVfWjJCyKBLtNWGVlQZHeW2R9CLgOPZcXu4Tku6R\ndDfwW0Dbxk4XXoJeZVclIfSTqtxjnmpuslMjvq5eHaTXroyjKD0ze9rM9jezH8fS3mNmrzKzXzez\nE2PRXiouvBh3Pjxnwh90nj/wqkhgUFRJ7r1Ge9vnTapstJe3XKMovaK48BLkHSjsVOfzySu9OlRz\nuy2HS687vNMiEFVlXXZNJzlDI9mpAcMZwtKrcJPjDp32eISHr1hchKpUb7tdY6/dczNaEV+/o76i\n12mV94gP+iyMvIx8hNdLJ4UzkfhnNqyByt1OR0sOYYE9Iz7YM/LKivy6Ge5ShkTTorojPriE9R/x\n2RidGPkILxmhuOyKM+zPsNtn4XaK+FrEI7JkdJZHZP2OGD3S68xIR3h3Pjzn+ffD/iN1yqeXiC9t\n2akWnZafahfplSU7b7MrxshGeHHZOc2mjOXkIT3qS5IUW1ltgXmuDR7ldWKkIzynf1ShTS9Jt8vJ\nw8QFCSA7ykpr9+sGj976z0gKz6O7wZKnuWCQUsyzwGiLduJL0k31t59y8w6MbEqp0ko6XtJ9kjZK\nOi9l/3slbZO0OrzOiO1bJGlDeC0qozztcNlVk0G3ofbasdGuyhsnKbRWlTRv1bQoXrVNp3CEF5Zb\n/jTRpN5NwF2Slqc8fexqMzs7cexU4AJgjGg9q1Xh2CeLlsupH2lPXes33UR7ceLSy4r+vIpaPcqI\n8BYAG83sATN7BrgKyLfqH7wNWGFmTwTJrQCOL6FMqXh0Vw+qHu0lyRv19Yu0yHPH4TuYc8XFQytT\nVSmjDe8g4JHY9ibgyJR875R0FPBD4BwzeyTj2NQWX0mLgcUAs2bNKqHYTpUZVrQH3T1HI06W9PK0\n/5XFMMVbB8oQXtpTNJLrxn8DuNLMdko6E1hGNHE1z7FRYrQC6lKIlnjvtpAe3dWTYfT2FhVfkk4S\nyiPEXkU254qLefC083s6tomUUaXdBBwS2z6YaM355zGz7Wa2M2x+AXhd3mPLwGXn9ELRqm5e4tXR\nZOdIGdVlr9rupgzh3QXMlTQnPGDjZGB5PIOkmbHNE4H14f2NwFslTQlPG3prSHOcytASX7/l59XR\n/lNYeGb2LHA2kajWA18zs7WSLpR0Ysj255LWSvoB8OfAe8OxTxA9aeiu8LowpDnOBKow/W9QUV/Z\neJQXUcrAYzO7Hrg+kfah2PvzgdSGBDO7DLisjHKkEV/nrgp/ME4xzpqxshIzN+LSK6utz+k/IzXT\noskP2hklhtGD245kxFdVAXoHRsMXD8ha2LMqfyhOMaqy+GiSulZ7R4HGCq/dKsZV/CNxeqfK4qsa\no96WN1JVWqfZVK2qC/Wp7o4KjRSeR3ejTRXF16IKAhzltrzGVmkdpw7/3AY1xs+JaGSEl0Ud/gCc\ncklGex79jTaNi/D8kYtOGnV8WFNa5FdWNDiqnReNE57j5KUO0oP0aq9XgXtjZIRXly+3M1jq/L2I\ni3AUBChpP0nXSLpX0npJb5A0VdKKsGL6ijAnP5NGCc+rs04vtKq7VR3Pl5c06bUTYQ2rtf8AfNvM\nDgdeTTR3/zzgJjObC9wUtjNpZKdFlRumnepT5+9PlvTq3gEi6SXAUexeeOQZ4BlJC4G3hGzLgFuA\nc7PO00jhOU4Z1Fl8SbKk188xeT/ZObkb0U6TNB7bXhoW/W1xKLAN+LKkVwOrgPcBM8xsC4CZbZF0\nQLuLNEZ4Xp11+kVVVmgpSsVXeHnczMba7N8LeC3wZ2Z2p6R/oEP1NY1GteFB+n/lOrfLONWg7u17\nSWrYybEJ2GRmd4bta4gEuLW1wHD4+Vi7kzROeI7TT5omvaPmbqxF54WZ/Qh4RNLLQ9IxwDqi1dVb\nz7NeBFzX7jyNEF6rOtukL6NTXZr2PatRtPdnwBWS7gbmAx8DLgGOk7SB6NnYl7Q7QWPa8OJ4ddbp\nN01p16sTZrYaSGvnOybvOUqJ8CQdL+k+SRslTWhIlPR+Sesk3S3pJkkvi+3bJWl1eC1PHus4jlMW\nhYUnaRLwaeAEYB5wiqR5iWz/BYyZ2a8TNTZ+Irbv52Y2P7xOpEuSvbP+X9cZFE2qPYzKKIcyIrwF\nwEYzeyAMBrwKWBjPYGY3m9nTYfMOoufPlkral69JX0jHcYpThvAOAh6JbW8KaVmcDtwQ254saVzS\nHZJOyjpI0uKQb3zbtm3FSuw4JeH/VOtFGZ0WSkmz1IzSHxA1Or45ljzLzDZLOhRYKekeM7t/wgmj\nUddLAcbGxlLP7zhVokkzNZpCGRHeJuCQ2PbBwOZkJknHAn8DnGhmO1vpZrY5/HyAaB7ca/Je2Nvv\nnCqQjPLquPYejEY7XhkR3l3AXElzgEeBk4FT4xkkvQb4PHC8mT0WS58CPG1mOyVNA97Inh0ajlML\nOknNo71qUDjCM7NngbOBG4mWa/mama2VdKGkVq/r3wO/DPxLYvjJEcC4pB8ANwOXmNm6XsrhXySn\nDtQl2msqpQw8NrPrgesTaR+KvT8247jvAq8qowyOUxd80PLwaMTUMsepG1WN9JrejjcywvP/qI7j\n1FZ4Tf9P5IwGVY30mkpthdeJ+BfJozunirS+oy69wdFY4TlOnWjaAqNVxYXnOBXGRVguLjzHqRB1\nnKFRJxovPG+/c5qAy68cGi+8Fv6FceqCf1f7Ry2F99DPHuqYxyM7x+mNJg/5qqXw8uLSc+qKL2jb\nHxotPMdxnDgjIzyP9pym4ENVemdkhAeR9Fx8Tp1JG7bi8svPSAnPcZqKSy8fLjzHcUYGF57jNASP\n8jrjwnMcZ2QoRXiSjpd0n6SNks5L2b+PpKvD/jslzY7tOz+k3yfpbWWUx3Gc5iFpkqT/kvTNsH25\npAfDc3JWS5rf6RyFn2khaRLwaeA4okc23iVpeeJhPKcDT5rZYZJOBj4O/L6keURPOXsFcCDwHUm/\nZma7ipbLcZzG8T6iB4W9JJb2l2Z2Td4TlBHhLQA2mtkDZvYMcBWwMJFnIbAsvL8GOEaSQvpVZrbT\nzB4ENobzOY7jPI+kg4HfBr5Y5DxlCO8g4JHY9qaQlponPNbxx8D+OY8FQNJiSeOSxnc8taNQgX0s\nntMU0h4CXlOmtf6+w2txYv+lwF8BzyXSL5J0t6QlkvbpdJEyHtOolDTLmSfPsVGi2VJgKcC0I6al\n5nEcpzq8YIeYfO/kvNkfN7OxtB2Sfgd4zMxWSXpLbNf5wI+AvYnccC5wYdsy5S1NGzYBh8S2DwY2\nZ+WRtBfwUuCJnMf2FY/2nKZR4ygvizcCJ0p6iKjJ7GhJ/2xmWyxiJ/BlcjSHlSG8u4C5kuZI2puo\nE2J5Is9yYFF4/y5gpZlZSD859OLOAeYC3y+hTB2Ji86l5zSNJknPzM43s4PNbDaRX1aa2R9ImgkQ\n+gNOAtZ0OlfhKq2ZPSvpbOBGYBJwmZmtlXQhMG5my4EvAV+VtJEosjs5HLtW0teAdcCzwJ8OsofW\nRec4teYKSdOJmsZWA2d2OqCMNjzM7Hrg+kTah2LvdwDvzjj2IuCiMspRhM9sPbpR/xUdp4mY2S3A\nLeF91xGLz7RwHGdkcOG1wau8jtMsXHiO44wMLrwYHtE5TrNx4TmOMzK48BJ4lOc4zcWFl4JLz3Ga\nSS2FN/vFs/t+jZb04j/9IUDOKLBswZeGXYS+UUvhDZqk5Fx6jlNPSplp4ThO9YjPHPJ/0hEe4TmO\nMzK48BxnBBnVeeMuvB7xKoJTV0ZVdlBj4VWhJ8ml59SNUZYdeKdFYeLSG/Uvk+NUndpGeFWk0zg9\njwidqpD1XaxCzamfuPD6QHLQcouzZqx06TnOEHHh9QmP9Jyq4U0uBYUnaaqkFZI2hJ9TUvLMl/Q9\nSWvD8yN/P7bvckkPSlodXvOLlKdOuPQcZ/AUjfDOA24ys7nATWE7ydPAH5rZK4DjgUsl7Rfb/5dm\nNj+8VhcsT+WIt+v5FDVnGPj3bDdFhbcQWBbeLyN6VNoemNkPzWxDeL8ZeAyYXvC60QVr1MCa9aXz\nL6MzCPJMM6vT31OvFBXeDDPbAhB+HtAus6QFRE8Jvz+WfFGo6i6RtE+bYxdLGpc0vm3btoLFrhYu\nPccZDB2FJ+k7ktakvBZ2c6Hw0NyvAn9kZs+F5POBw4HfAKYC52Ydb2ZLzWzMzMamTy8lQKwULj2n\nX3hnxW46Cs/MjjWzV6a8rgO2xp7+PZOoujoBSS8BvgX8rZndETv3FovYCXwZWNDtDTQpDHfpOYPg\nrBkrJ0iwSX9H7ShapV0OLArvFwHXJTNI2hv4OvAVM/uXxL6WLEXU/remYHkcx3EyKSq8S4DjJG0A\njgvbSBqT9MWQ5/eAo4D3pgw/uULSPcA9wDTgowXLU3s8ynOc/lFoLq2ZbQeOSUkfB84I7/8Z+OeM\n40v561624Ess+v7pZZyqEvj8XKffjOo/1sYsHnDbhsP22D5q7sYhlaRcPrP1aJee01fq0H4naTJw\nG7APkbeuMbMLJM0BriLq9PxP4D1m9kzWeXxqWQ1otyjBqP6ndkaOncDRZvZqYD5wvKTXAx8HloTJ\nD08Cbat6jYnwkty24bDGRHktOknPI0GnqZiZAT8Nm78UXgYcDZwa0pcBHwY+m3WexkR4D552/rCL\nMHQ82nPyEP+eVKg6O601sSC8FiczSJokaTXR8LcVRBMYnjKzZ0OWTcBB7S7S2AgP9ozymhjxOU6V\nmfRz2H/drrzZHzezsXYZzGwXMD/Mxf86cERatnbnaLTwYM/OjPj7JsrPq7RON9T1+2JmT0m6BXg9\nsJ+kvUKUdzCwud2xjanSQnfV2mSv7m0bDnv+5ThNJV6dPfJlDw6xJN0haXprlSVJLwSOBdYDNwPv\nCtlSJz/EaZTwuiVLcq3tusnP2/CcBjMTuFnS3cBdwAoz+ybR/Pv3S9oI7A+0bZRsXJX2wdPOZ84V\nFxc+TyfZtfZXrWrsPbZOFvHvRp2iOwAzuxt4TUr6A3QxB3+kI7w8dKrmVjUK9GjPcSbSuAgPyovy\n4nQrvSpEfj5Lw4nT+i58ZuvRHPmyIRdmSHiE1yeq0gHS6dGRzuhRobF3A6exwqvKQOQqic9x6tZ2\nVzaNrNJWkSqMAfQODWfUaWyEV2WGPebPoz1nVGm08KpSrW3HsMTnbXujx6hXZ6HhwoN6SA+GN7yl\nJb60Z+a6EJ2m4W14FWLYi5imCS6Z5u1/9cSju4hCEZ6kqZJWSNoQfk7JyLcr9jyL5bH0OZLuDMdf\nHR744wSq0sPrOE2haJX2POCmsNroTWE7jZ+b2fzwOjGW3tVqpb3SS7V28r2Tn39lpXd6lYVLzymC\nR3e7KSq8hUSrjBJ+npT3wPBoxqOBa3o5vgzySqlXiZUpwKpEe96u59SZom14M8xsC0QP1ZZ0QEa+\nyZLGgWeBS8zs34hWNsi9WmlYAXUxwKxZs7ouaHy6WVw+ZUZinci61o7Dd+Q+RxUWLfDxfNWh0+/C\no7s96RjhSfqOpDUpr4VdXGdWWM30VOBSSb8KKCVf5mqlZrbUzMbMbGz69OldXHo3D552/kAFl5de\nokCP9pz455/Wq+6ym0jHCM/Mjs3aJ2mrpJkhuptJtNZ82jk2h58PhJVKXwNcS5erlVaR+BLW2+dN\nKny+lvTyRH0+e8NJ4gtGtKdoG95yolVGIWO1UUlTJO0T3k8D3gisC08h6mq10jJY/5FzSjnP/ut2\nTVivv5WWfPVCLxHfMKM+H7c3WDp91h7dpVO0De8S4GuSTgf+G3g3gKQx4EwzO4PoQRufl/QckWAv\nMbN14fhzgaskfRT4LzqsVloW6z9yDkd8cElPx/YisG6OSUaJ3UR8MPw2vrzS8yikd9p9xnVc3HOQ\nFBKemW0HjklJHwfOCO+/C7wq4/iuVisdNp3Ete+qRwH4yevaPimu4zXSqsZ1E18nvOrVG51kN8pr\n3eWh8VPLsuimapunatqSXet9fLtb2lWH61bVdcojT/Q8ymvd5cGnlrUhT1W0ndg6RXztjm0dkyxD\nK/rrJeKrarTndMbbR8thZCM8yI7y8kZ0eaO4eN5uj0nmT5arm4jPo716kkd23naXj5GP8OIdGL32\nqOalW+klj9tNFP3F2/p6Gc4y7IjPh7R0xmVXLiMd4bXoZvhI0fa5MmhdP63c3UZ8TnXJKzsnPy48\n4PZrP5Ar37BFFycu3qzOjTx4NbeapC3L1Xol83l0lx8XXqCd9KoQ1WXRKludoz1vkN+TvGsQnjVj\npffKdokLL0ZSelUWXZIyqrke7Q2fTvKP7/fIrntceAnyVm+rSJnVXGfwpMkuK7rzyK43XHgp3H7t\nB2oT2aXRqZqbh0FLr+lzcbuJ3FqktdeNKpIuk/SYpDWxtA9LejS2mvrbO51n5IelNJnd0t5zGMvk\neydXdvhKE4eqlC27EY3uLgc+BXwlkb7EzD6Z9yQe4WVww0O9LS5QRdIivqp3aGQ9Ta1udBK4yy4f\nZnYb8ETR87jw2tAk6cWJV3N9+Ep/iMvaZddXzpZ0d6jypj5ELI6iZenqxdjYmI2Pjw/seifMLmcN\nvSoQn9cbn6nRzTLzw5ih0amK202Dfz/JO6Sk6rKTtCqsUt4zL91nhv3mzFNz5f32w5c+DDweS1pq\nZksTZZoNfNPMXhm2Z4RjDPgIMNPM/rjddbwNLwc3PLSkMdLbd9WjqYsZVH1qWtWqtsnorRvh5rmX\nEYzsHu9WsGa2tfVe0heAb3Y6xoWXkyZJr0Xa2nt5OzSgOnNyixAXS56IMCmrLHl1K7t4/hGUXU+0\nHi8RNt8BrGmXH1x4XdEU6cUjTpObAAAIZ0lEQVSjvKLSqzJZi4xmSafdoqSdRJVHnO0iu7R9Lrvd\nSLoSeAswTdIm4ALgLZLmE1VpHwL+T6fzuPC6pKnSg4mrrzRBenF6rRZ3qq4WlV0aLrs9MbNTUpK7\n/pAK9dJKmipphaQN4eeEXhJJvxUbGLha0g5JJ4V9l0t6MLZvfpHyDIqm9N4mp87VZZByN3Q7vKXT\now+TE/i7bY/Lg8uufxQdlnIecJOZzQVuCtt7YGY3m9l8M5sPHA08Dfx7LMtftvab2eqC5RkYTZEe\nMEF6oz5sJUuQeWY+dBJmJ1x2/aWo8BYCy8L7ZcBJHfK/C7jBzJ4ueN1K0FTpFaFJ4mtHp/a4XqrP\nLrv+U1R4M1q9JOHnAR3ynwxcmUi7KAwcXNJ6fm2daJL04iSjvG4eHATVrub2Sj9nf7jsBkNH4Un6\njqQ1Ka+F3VxI0kyixzXeGEs+Hzgc+A1gKtFzarOOXyxpXNL4tm3burl032mC9PI8WtKl1x9cdoOj\no/DM7Fgze2XK6zpgaxBZS2iPtTnV7wFfN7NfxM69xSJ2Al+mzTNqzWypmY2Z2dj06dPz3t/AaIL0\n8tBLtOdk47IbLEWrtMuBReH9IuC6NnlPIVGdjclSRO1/HQcOVpkbHlpSW/H1azksj/KycdkNnqLC\nuwQ4TtIG4LiwjaQxSV9sZQpz4A4Bbk0cf4Wke4B7gGnARwuWpxLUVXpJ2j3YqNuHgTu7WbbgSy67\nIVFo4LGZbQeOSUkfB86IbT9Ea1G2PfNVa4JkiTRlgHLaTAynd1x0w8WXh+ojdYv0sqq1WZGeR3nd\n4bIbPi68PlO3dj2XXvl4FbY6uPAGRJ2k109GTXouumrhwhsgdZFeP6M8GB3pueyqh6+WMmBa0qt6\nh0ZLenkGJEP3q6vctuGwWq+j1w4XXXVx4Q2JOokvLr0ye23jkV4T5Oeiqz5epR0ydajmJqu47cbn\njSouu3rgEV4FqEu0F6fs8Xl1reK66OqFC69CVFl8WQ//idPEVZKzcNHVExdeBamy+DrRdOm56OqN\nC6/CVEl8eXtri8iu6p0YLrv6450WNaAKHRv9Wk2lDvhMiebgEV5NiEuvChFfFmVVaYf9zFsXXDNx\n4dWQYVR181Zpodx2vEFXc110zcaFV2OqHPXVqfPCJTc6uPAaQpU6OPpF2WP1XHSjhwuvYVQp6utH\nlFdUei650caF12DKkF9W212dVkF2yTktXHgjQnJoSx4BdtNRMUg6dWS44JwsCglP0ruBDwNHAAvC\nsyzS8h0P/AMwCfiimbUe9jMHuIrombT/CbzHzJ4pUiYnH50EWJbs+t15cduGw3jwtPP7dn6nOmR5\npBuKRnhrgN8FPp+VQdIk4NNETzXbBNwlabmZrQM+Diwxs6skfQ44HfhswTI5PZA2uPlN7/xkKecu\nW3ouuNGjg0dyU/SpZetDYdplWwBsNLMHQt6rgIWS1gNHA6eGfMuIokUXXkW4/doPdMxzxAc7zwLJ\nKzsXmdOGVI8AgxNeTg4CHoltbwKOBPYHnjKzZ2PpmfUoSYuBxWFzp6RaP7Q7g2nA48MuRJ/oeG/6\ng78eUFFKpam/s5cXPcH/PPPYjd9++NJpObNPlhRvEltqZktj21ke6YqOwpP0HeBXUnb9jZldl+Ma\naeGftUlPJdz80lCmcTMby3HtWtHU+4Lm3luT76voOczs+DLKEujKF1l0FJ6ZHdvtSRNsAg6JbR8M\nbCb6r7ifpL1ClNdKdxzHSZLlka4YxGopdwFzJc2RtDdwMrDczAy4GXhXyLcIyBMxOo4zeqR6pNuT\nFBKepHdI2gS8AfiWpBtD+oGSrgcI0dvZwI3AeuBrZrY2nOJc4P2SNhK16eUdQLW0c5Za0tT7gube\nm9/XAOjgkdwoCrQcx3Gajy8A6jjOyODCcxxnZKiF8CS9W9JaSc9JyhwCIOl4SfdJ2ijpvEGWsRck\nTZW0QtKG8HNKRr5dklaHV9cNtYOi0+cvaR9JV4f9d0qaPfhS9kaOe3uvpG2x39MZwyhnN0i6TNJj\nWWNaFfGP4Z7vlvTaQZexdMys8i+iubovB24BxjLyTALuBw4F9gZ+AMwbdtk73NcngPPC+/OAj2fk\n++mwy5rjXjp+/sBZwOfC+5OBq4dd7hLv7b3Ap4Zd1i7v6yjgtcCajP1vB24gGgP3euDOYZe56KsW\nEZ6ZrTez+zpke37qiUULELSmnlSZhURT6gg/TxpiWYqS5/OP3+81wDHqMC+xItTxu9URM7sNeKJN\nloXAVyziDqJxszMHU7r+UAvh5SRt6kk11zfazQwz2wIQfh6QkW+ypHFJd0iqqhTzfP7P57FomMGP\niYYjVZ283613hqrfNZIOSdlfN+r4N9WWyqyH18cpbEOl3X11cZpZZrZZ0qHASkn3mNn95ZSwNPJ8\n/pX8HeUgT7m/AVxpZjslnUkUyR7d95L1l7r+vjKpjPCsf1PYhkq7+5K0VdJMM9sSqgqPZZxjc/j5\ngKRbgNcQtSlViTyffyvPJkl7AS+lfZWqKnS8NzPbHtv8AtHSZ3Wnkn9TRWhSlbaUqScDZjnRlDrI\nmFonaYqkfcL7acAb6XJJnAGR5/OP3++7gJUWWscrTsd7S7RtnUg0G6DuLAf+MPTWvh74casJprYM\nu9ckZ2/SO4j+2+wEtgI3hvQDgesTvUo/JIp+/mbY5c5xX/sDNwEbws+pIX2MaEVXgN8E7iHqGbwH\nOH3Y5W5zPxM+f+BC4MTwfjLwL8BG4PvAocMuc4n3djGwNvyebgYOH3aZc9zTlcAW4Bfh7+t04Ezg\nzLBfRItu3h++e6kjJOr08qlljuOMDE2q0jqO47TFhec4zsjgwnMcZ2Rw4TmOMzK48BzHGRlceI7j\njAwuPMdxRob/D1E3OGXfex++AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x103f289e8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import scipy.interpolate\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "x = np.linspace(-1,1, num=500)\n",
    "y = np.linspace(-1,1, num=500)\n",
    "grid_x, grid_y = np.meshgrid(x,y)\n",
    "\n",
    "X = res[:,0]\n",
    "Y = res[:,1]\n",
    "disp_x = scipy.interpolate.griddata((X, Y), res[:,2], (grid_x, grid_y), method='linear')\n",
    "plt.figure()\n",
    "plt.gca().set_aspect('equal')\n",
    "plt.contourf(grid_x, grid_y, disp_x)#,levels=np.linspace(75,180,20))\n",
    "plt.colorbar();\n",
    "plt.title(\"Breakout width\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, the safest recommentation would be to drill horizontal in the direction of $S_{Hmax}$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Extra test cases"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 43.08493649,  -4.82962913,   0.08493649],\n",
       "       [ -4.82962913,  46.83012702,  -4.82962913],\n",
       "       [  0.08493649,  -4.82962913,  43.08493649]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S = np.array([[50, -5, 0],[-5, 40, 0],[0, 0, 43]]);\n",
    "compute_Sb(S,(0,0,0),(30,45))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(44.464101615137757, 97.320508075688764, 0.0, 0.0)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compute_wellbore_stress(S, 0.2, 30, 0.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "97.320508075688764"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compute_max_tangent_stress(S, 0.2, 30, 0.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "44.464101615137764"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compute_min_tangent_stress(S, 0.2, 30, 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
}