scarp.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('wet angle of repose: ', 0.46630765815499858, 25.0)\n",
      "('beach slope: ', 0.02, 1.1457628381751035)\n",
      "('failure angle: ', 1.0076445808463947, 45.218166156499294)\n",
      "('dune slope: ', 0.05, 2.8624052261117479)\n",
      "('Eroded volume: ', 0.18084710838307211, 0.020152891616927897, 0.0010000000000000002)\n",
      "('A,B,C', array([  23.85423716,  112.07335952,   44.07240332]))\n",
      "('b, c = ', 1.0915914854264472, 0.81932879156631067)\n",
      "('xs, xl', 0.81916497495134399, 0.17015288855354893)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1f76208>]"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHbBJREFUeJzt3Xlw3PWZ5/H3I8myJV+SLNkyOpDxgQ9BgCjGgQDmCMgs\nA0kIGcg9y8S7mZCdSTKzIckuyTCVqkyytTOzGzKsK6Gyyc5CYeXASXxwhIQlxAkmEJBtDoXDkm1Z\ntiVf6LCOZ//oltySW+qW1K3u30+fV5XLfXy7+4FWf/zTt5/v92fujoiIhEtOpgsQEZHUU7iLiISQ\nwl1EJIQU7iIiIaRwFxEJIYW7iEgIKdxFREJI4S4iEkIKdxGREMrL1AuXlpZ6TU1Npl5eRCSQnnvu\nuSPuXpZoXMbCvaamhl27dmXq5UVEAsnM3kpmnKZlRERCSOEuIhJCCncRkRBSuIuIhJDCXUQkhBKG\nu5k9YGZtZtaYYNy7zKzfzD6YuvJERGQikjly/z5QP9YAM8sF/hHYkYKaRERkkhKGu7s/BbQnGPZZ\n4EdAWyqKEhEJo7eOvs39v/4Tz/zpSNpfa9KLmMysAng/cA3wrgRjNwIbAaqrqyf70iIiWa+p7STb\nXmpla2Mrew+eAODT65dy2dLStL5uKlao/jPwRXfvN7MxB7r7JmATQF1dnc7MLSKh4+7sOXiC7Y2t\nbGtspantFACXVBfxlRtXUV9bTlVJYdrrSEW41wEPRYO9FLjRzPrc/acpeG4Rkazn7rzQfIztja1s\n393KW0c7yTFYu6SEj61bww1ryimfP2tKa5p0uLv7ksHLZvZ94OcKdhEJu/4B57m3OtjWeJAdja0c\nON5NXo7x7qUL+A9XLuX6NYsonTMzY/UlDHczexBYD5SaWQvwVWAGgLvfn9bqRESySF//ADtfb48E\n+u5DHDnVQ35eDlcuL+Xz15/PdasWUlSYn+kygSTC3d3vSPbJ3P2Tk6pGRCTL9PT180zTUbY1HuSx\nPYfo6OylYEYuV68so752MdesXMicmRnbYHdU2VeRiEiGdff286tXDrO98SBP7G3jZE8fc2bmce2q\nhWyoLeeqFQspyM/NdJljUriLiACnevp48uU2tje28suX2+jq7aeocAb1teVsuKCcy5eVMjMvuwM9\nlsJdRKat4529PL73ENsaW3nqtcOc7hugdE4+77+kghtrF3PpeSXMyA3mFlwKdxGZVo6e6uGxPZFA\n/03TEfoGnMXzZ/HhtdVsqC2nrqaE3Jyx1+wEgcJdRELv0IluduxuZdtLrfzujaMMOFSXFHLne5ZQ\nX1vOOyqLyAlBoMdSuItIKLV0dA6tEv3Dvg7cYWnZbP5q/TLqa8tZc848Eq2qDzKFu4iExhtH3mZb\n40G2N7byYstxAFYtnsfnrlvBhtpyli+am+EKp47CXUQCy915re0UW1+KBPrLrScBeEflfO7esJL6\nNeXUlM7OcJWZoXAXkUBxd3YfOMG2xoNsa2zl9cNvYwZ15xbzX29aTX1tORVFBZkuM+MU7iKS9QYG\nnOebj7E9GugtHV3k5hjrzivhLy5fwg2rF7Fw3tRuzJXtFO4ikpX6B5xn32yP7LTY2ErriW5m5BqX\nLyvls9cs472ryymZnR37uGQjhbuIZI3e/gF++6ejbGts5bE9rRw5dZqZeTlctaKML15wPtesXMT8\nghmZLjMQFO4iklHdvf08/doRtjW28vjeQxzv6qUwP5erVy7kxtrFrD+/jNlZuDFXttP/MRGZcp2n\n+/j1K4fZFt3H5VRPH3Nn5fHeVYuory3nyhVlzJoRnH1cspHCXUSmxMnuXn75chvbXmrlV6+20d07\nQMnsfG66cDH1teVctrSU/Lxg7uOSjRTuIpI2xzpP89ieQ2xvbOX/vXaE0/0DLJw7k9veWcWG2nLW\nLikhL6Abc2U7hbuIpNThkz08uifS4fLbPx2lb8CpKCrgY+8+lw215VxSXRy6fVyykcJdRCat9Xg3\n2xsPsrWxlV1vtjPgULOgkE9deR4basu5oGJ+qPdxyUYKdxGZkOb2zqFVos/vOwbAikVzuOua5Wyo\nLWdl+VwFegYlc4LsB4CbgDZ3r41z/0eAL0avngI+7e5/TGmVIpIVmtpODa0S3X3gBAC1FfP4uxvO\np762nKVlczJcoQxK5sj9+8C3gR+Mcv8bwFXu3mFmG4BNwKWpKU9EMsndebn1JNsaW9neeJBXD50C\n4OLqIr5840o21C6mqqQww1VKPAnD3d2fMrOaMe5/JubqTqBy8mWJSCY1tZ2k4bn9bG88yJtHO8kx\neFdNCV/7s9XcUFvO4vnamCvbpXrO/U5g22h3mtlGYCNAdXV1il9aRFLlo9/9PUdO9fDupQvYeOVS\nbsr5DfPan4HLv57p0iRJKWswNbOriYT7F0cb4+6b3L3O3evKyspS9dIikmJdvf185NJqfnjnpXz4\n0mrmte2CPz6Y6bJkHFJy5G5mFwLfBTa4+9FUPKeIZJY6XYJt0kfuZlYN/Bj4mLu/OvmSRERkspJp\nhXwQWA+UmlkL8FVgBoC73w/cAywAvhP9l77P3evSVbCIiCSWTLfMHQnu/0vgL1NWkYiITJp27BGR\nxDZ8C/62KdNVyDho+wERSSxHx4FBo3dMRBJ7qQF2fCXTVcg4KNxFJLF9v1Wfe8Ao3EVEQkjhLiIS\nQgp3EZEQUriLiISQwl1EElOfe+Coz11EElOfe+DoHRORxNTnHjgKdxFJTH3ugaNwFxEJIYW7iEgI\nKdxFREJI4S4iEkIKdxFJTH3ugaM+dxFJTH3ugaN3TEQSU5974CQMdzN7wMzazKxxlPvNzP6HmTWZ\n2YtmdknqyxSRjFKfe+Akc+T+faB+jPs3AMujfzYC/zr5skREZDIShru7PwW0jzHkFuAHHrETKDKz\nxakqUERExi8Vc+4VQHPM9ZbobSIikiGpCHeLc5vHHWi20cx2mdmuw4cPp+ClRUQknlSEewtQFXO9\nEjgQb6C7b3L3OnevKysrS8FLi8iUUJ974KQi3LcAH492zawDjrv7wRQ8r4hki5wc9boHTMJFTGb2\nILAeKDWzFuCrwAwAd78f2ArcCDQBncBfpKtYEcmQlxrgwPNww9czXYkkKWG4u/sdCe534DMpq0hE\nss++38LunyjcA0S/Z4mIhJDCXUQkhBTuIiIhpHAXEQkhhbuIJKY+98DRfu4ikph63ANH75iIJKb9\n3ANH4S4iiWk/98BRuIuIhJDCXUQkhBTuIiIhpHAXEQkhhbuIJKY+98BRn7uIJKY+98DROyYiianP\nPXAU7iKSmPrcA0fhLiISQgp3EZEQUriLiIRQUuFuZvVm9oqZNZnZ3XHurzazJ83seTN70cxuTH2p\nIiKSrIThbma5wH3ABmA1cIeZrR4x7L8AD7v7xcDtwHdSXaiIZJD63AMnmSP3tUCTu7/u7qeBh4Bb\nRoxxYF708nzgQOpKFJGMy8lRr3vAJPNuVQDNMddborfF+hrwUTNrAbYCn01JdSKSHdTnHjjJhLvF\nuc1HXL8D+L67VwI3Aj80s7Oe28w2mtkuM9t1+PDh8VcrIpmhPvfASSbcW4CqmOuVnD3tcifwMIC7\n/xaYBZSOfCJ33+Tude5eV1ZWNrGKRUQkoWTC/VlguZktMbN8Il+YbhkxZh9wLYCZrSIS7jo0FxHJ\nkITh7u59wF3ADmAvka6Y3WZ2r5ndHB32BeBTZvZH4EHgk+4+cupGRESmSFK7Qrr7ViJflMbedk/M\n5T3A5aktTUREJkq9TSKSmPrcA0f7uYtIYupxDxy9YyKSmPrcA0fhLiKJqc89cBTuIiIhpHAXEQkh\nhbuISAgp3EVEQkjhLiKJqc89cNTnLiKJqc89cPSOiUhi6nMPHIW7iCSmPvfAUbiLiISQwl1EJIQU\n7iIiIaRwF5HE+gcyXYGMk8JdRMb28svwyf8NF92f6UpkHBTuIjK2hgbo6IA1azJdiYyDwl1ExrZ5\nM1x8PjR+J9OVyDgo3EVkdK++Ci++CHWL1eceMEmFu5nVm9krZtZkZnePMuZDZrbHzHab2f9NbZki\nkhENDZG/L1ue2Tpk3BLuLWNmucB9wHuBFuBZM9vi7ntixiwHvgRc7u4dZrYwXQWLyBTavBnWrYOy\nudCW6WJkPJI5cl8LNLn76+5+GngIuGXEmE8B97l7B4C768dAJOBKDu6DF16A227LdCkyAcmEewXQ\nHHO9JXpbrBXACjP7jZntNLP6eE9kZhvNbJeZ7Tp8+PDEKhaRKbH6N49GLnzwg5ktRCYkmXC3OLf5\niOt5wHJgPXAH8F0zKzrrQe6b3L3O3evKysrGW6uITKHVzzwKa9dCdbX2cw+gZMK9BaiKuV4JHIgz\n5hF373X3N4BXiIS9iARQRftBzvnTnjNTMjk52tM9YJJ5t54FlpvZEjPLB24HtowY81PgagAzKyUy\nTfN6KgsVkalz3d6nIxduvTXyt/ZzD5yE4e7ufcBdwA5gL/Cwu+82s3vN7ObosB3AUTPbAzwJ/J27\nH01X0SKSXtfvfor9y9bAkiWRG7Sfe+AkdZo9d98KbB1x2z0xlx34fPSPiATZW29Re+BVHvv4587q\nnJDg0CSaiAwXXbi057L3ZrgQmQyFu4gM19DA3vJldCyuznQlMgkKdxE5o7kZdu7k0TVXZLoSmSSF\nu4ic8aMfAfD46hHhrj73wFG4i8gZmzfDO97BvgUjvkpVn3vg6N0SkYj9++GZZ+JvN6A+98BRuItI\nRHRKJu5GYepzDxyFu4hENDTABRfA+ednuhJJAYW7iMDBg/D009oBMkQU7iICP/4xuGvv9hBRuItI\npEtm9WpYtSrTlUiKKNxFprtDh+Cpp8Y+alefe+AktXGYiITY4JTMWPPt6nEPHL1jItPd5s2wciWs\nWTP6GPW5B47CXWQ6a2uDX/86ctRu8c6oGaU+98BRuItMZz/5CQwMqEsmhBTuItNZQwMsXx5ZvCSh\nonAXma6OHIEnn4wctY81JSOBpHAXma5++lPo79eq1JBKKtzNrN7MXjGzJjO7e4xxHzQzN7O61JUo\nImmxeTMsXQoXXZR4rPrcAydhuJtZLnAfsAFYDdxhZqvjjJsL/Cfgd6kuUkRS7OhReOKJ5KdktJ97\n4CTzbq0Fmtz9dXc/DTwE3BJn3D8A3wS6U1ifiKTDI4+Mb0pGfe6Bk0y4VwDNMddborcNMbOLgSp3\n//lYT2RmG81sl5ntOnz48LiLFZEUaWiAJUvgkkuSG//ao7DnkfTWJCmVTLjH+53Nh+40ywH+CfhC\noidy903uXufudWVlZclXKSKp09EBjz+eeOFSrCOvwYKl6a1LUiqZcG8BqmKuVwIHYq7PBWqBX5nZ\nm8A6YIu+VBXJUlu2QG9v8guX3OFoEyxYnt66JKWSCfdngeVmtsTM8oHbgS2Dd7r7cXcvdfcad68B\ndgI3u/uutFQsIpOzeTOcey7UJXn8daoNek5AqcI9SBKGu7v3AXcBO4C9wMPuvtvM7jWzm9NdoIik\n0PHj8Oij45uSOdECMwphwbL01iYpldSWv+6+Fdg64rZ7Rhm7fvJliUhaDE7JjGfhUsU74Uv7ifmq\nTQJA+7mLTCcNDVBVBZdeOr7Hqcc9cPSOiUwXJ07Ajh1w663j20vmsXvgV99IX12SFgp3keni5z+H\nnp7xb++792fQtjc9NUnaKNxFpovNm6GiAtatS/4xfT3Q8aY6ZQJI4S4yHZw8Cdu2RaZkxjN/3v4G\n+ACUrkhfbZIWCneR6eAXv5jYlMzR1yJ/qw0ycBTuItNBQwMsXgyXXTa+x/V2wZxyhXsAqRVSJOze\nfhu2boU77xx/S+OFH4r8kcDRkbtI2P3iF9DVpTMuTTMKd5Gwa2iARYvgPe8Z3+Pc4XvXw64H0lOX\npJXCXSTMOjsjR+4f+ADk5o7roYV9x6D5d5F5dwkchbtImG3bFgn48XbJAKU9+yIXtNVvICncRcJs\n82YoK4Mrrhj3Qxd0R8NdC5gCSeEuElZdXZEtBz7wAcgbf2Ncafc+yJ0JRdVpKE7STeEuElbbt0fa\nICfYJdOVNxeWXQs545url+ygcBcJq4YGWLAA1q+f0MOfLv843PFgamuSKaNwFwmj7m742c/g/e+f\n0JSMBJ/CXSSMHn00slnYBLpkAM71A3zupfdB0xMpLkymisJdJIw2b4aSErj66gk9vIb9zO9tg1nz\nU1yYTJWkwt3M6s3sFTNrMrO749z/eTPbY2YvmtkTZnZu6ksVkaT09ETOlfq+98GMGRN6ihoORC5o\nw7DAShjuZpYL3AdsAFYDd5jZ6hHDngfq3P1CoAH4ZqoLFZEkPfZY5JR6E5ySgUi4n8orhoKiFBYm\nUymZI/e1QJO7v+7up4GHgFtiB7j7k+7eGb26E6hMbZkikrSGBigqgmuumfBT1HCAo7PU3x5kyXyN\nXgE0x1xvAcY6dfqdwLbJFCUiE3T6NDzySGRKJj9/1GF9/QMcPN7NvvZOmts7ae7opLm9i33tnbR0\ndPJXfTXUzF+C5leDK5lwj3eadI870OyjQB1w1Sj3bwQ2AlRX66hAJOWeeAKOHcNvvZUjJ3uGwrq5\n/Ux4N3d0cvB4N/0DZz7GuTlGRVEBVSUFXLdqEV0lX2flxRUZ/A+RyUom3FuAqpjrlTD4bcsZZnYd\n8BXgKnfvifdE7r4J2ARQV1cX9x8IEUnOye5emtu7okfdnbR0dHHtt/6Vi2fN5rLfDHDi6ceHjS+b\nO5Oq4gLeeW4xVcWFVJcUUllSQFVxIYvnzyIvV81zYZJMuD8LLDezJcB+4Hbgw7EDzOxi4H8B9e7e\nlvIqRaah030D7D/WNWza5MzlTjo6e4eNnzMzjwtmFLLzyj/jtsuXUV1SSFU0vCuLCynI1zYC00nC\ncHf3PjO7C9gB5AIPuPtuM7sX2OXuW4BvAXOAzWYGsM/db05j3SKBNzDgtJ3sGQrr2GmTlvZODp7o\nxmN+v52Ra1QWF1JZXEDtBYuHjr4HA7yocAZmNwBwXYb+myR7JLUu2d23AltH3HZPzGX9LInEcbyz\n90x4d3RGv8CMTKW0dHRxum9g2PjyebOoKilg3XkLqCyJhndxAVUlhSyaN4vcnHhfgYmcTZtOiExC\nd28/LR1dQ0fbseHd3N7Jie6+YePnF8ygqqSA8xfN5bpVi6iKCe+KogJmzdDUiaSGwl1kDP0DTuuJ\n7ui0yeAR+Jm570MnhvcOzMzLoTIa1pdUFw9Nm1QWF1JVUsj8gomtGBUZL4W7TGvuTkdn71nTJoPt\ng/uPddHbf2biO8dg8fwCKosLuGJ52bA576qSQsrmzCRHUyeSBRTuEnqdp/siUydxpk2a2zt5+3T/\nsPEls/OpKi5gTcV8NkS/uBwM8HOKCsjPU8ugZD+FuwTe4GrLofDuGB7gR06dHja+YEbuUFivO2/B\nsHnvqpJC5szUx0KCTz/FkvXcnSOnTg872h4M733t8VdbnlM0i6riwqEvLQfnwatLClkwO59oy65I\naCncJSuc6umLmTbpHDaN0tLRRVfv8KmT0jkzqSo5s9oydt5bqy1FFO4yRU73DXDgWFfcaZPRVltW\nFhewpHQ2V64oG5o2qS7RakuRZCjcJSUGBpzDp3rO7DIYM23S0t5J64luBkastoxsVFU4tNqyqqQg\numhncLWlpk5EJkrhLkk73tUb0+s9Yrl8nNWWi+bNpLqkcGi1ZVVxNLy12lIk7RTuMqS7t5/90amT\nlpjFOoNH4yNXW86blUf1gsIzqy1jOk602lIksxTu00j/gHPoROwJGrrOLJmPs9oyPy9nKLAvqS4e\nmjbRakuR7KdwDxF351hn76hfWo5cbWkGi+fNoqqkkCuWlw2f99ZqS5FAU7gHTNfp/uH93h1dw9oH\nT/UMnzopLpxBdUkhayrmU1+7eNiXllptKRJeCvcsE7vacuSXls3tXRw5NXzqZNaMnKGwXnfeAipj\nvrTUakuR6Uuf/Ck2crVlS0cX+45Gw7ujkwPHRl9tee3KhZHFOoPhXVxI6RytthSRsync02BwteWw\n7WFjjr5HW215cVUxN7/jzLSJVluKyEQp3CdgcLXlWb3e0TBvf3v4RlWDqy3PXTCb9ywrozrm6Luy\nuIDCfL0NIpJaSpU4BldbDu3xfXT4GeYPHu8adbXlDefMP2uP72KtthSRKZZUuJtZPfAvRE6Q/V13\n/8aI+2cCPwDeCRwF/tzd30xtqak1uNqypSPOHt+jrLasKi5k7ZKSs7aILddqSxHJMgnD3cxygfuA\n9wItwLNmtsXd98QMuxPocPdlZnY78I/An6ej4GQNrraMd2q0fUfjr7asKilk+cK5XBtdbVkZnfuu\nLNZqSxEJlmSO3NcCTe7+OoCZPQTcAsSG+y3A16KXG4Bvm5m5u5Mmg6stz5ygoSs65x25Hm+1ZWVx\nZKrkoqqiYV9aVhUXMr9Qqy1FJDySCfcKoDnmegtw6Whj3L3PzI4DC4AjqSgy1pMvt/H3P9s96mrL\nypLC6JeWhWfaBosLWThXqy1FZPpIJtzjJeLII/JkxmBmG4GNANXV1Um89NmKZ+ez5pwzqy0Hj77P\nKZrFzDxNnYiIQHLh3gJUxVyvBA6MMqbFzPKA+UD7yCdy903AJoC6uroJTdlcVFXEfR+5ZCIPFRGZ\nNpJZHfMssNzMlphZPnA7sGXEmC3AJ6KXPwj8Mp3z7SIiMraER+7ROfS7gB1EWiEfcPfdZnYvsMvd\ntwDfA35oZk1EjthvT2fRIiIytqT63N19K7B1xG33xFzuBm5LbWkiIjJR2rRERCSEFO4iIiGkcBcR\nCSGFu4hICCncRURCyDLVjm5mh4G3JvjwUtKwtUEaBaneINUKwao3SLVCsOoNUq0wuXrPdfeyRIMy\nFu6TYWa73L0u03UkK0j1BqlWCFa9QaoVglVvkGqFqalX0zIiIiGkcBcRCaGghvumTBcwTkGqN0i1\nQrDqDVKtEKx6g1QrTEG9gZxzFxGRsQX1yF1ERMYQqHA3s9vMbLeZDZhZ3Yj7vmRmTWb2ipndkKka\nR2NmF5nZTjN7wcx2mdnaTNc0FjP7bPT/5W4z+2am60mGmf2tmbmZlWa6ltGY2bfM7GUze9HMfmJm\nRZmuaSQzq4++901mdnem6xmLmVWZ2ZNmtjf6s/rXma4pETPLNbPnzezn6XydQIU70Ah8AHgq9kYz\nW01km+E1QD3wneiJvbPJN4G/d/eLgHui17OSmV1N5Ly4F7r7GuC/ZbikhMysishJ3PdlupYEHgNq\n3f1C4FXgSxmuZ5jo5+Y+YAOwGrgj+vnKVn3AF9x9FbAO+EyW1wvw18DedL9IoMLd3fe6+ytx7roF\neMjde9z9DaCJyIm9s4kD86KX53P22ayyyaeBb7h7D4C7t2W4nmT8E/CfiXN6x2zi7o+6e1/06k4i\nZzbLJmuBJnd/3d1PAw8R+XxlJXc/6O5/iF4+SSQ0KzJb1ejMrBL4d8B30/1agQr3McQ7iXe2vcF/\nA3zLzJqJHAln1RHbCCuAK8zsd2b2azN7V6YLGouZ3Qzsd/c/ZrqWcfr3wLZMFzFCED5LcZlZDXAx\n8LvMVjKmfyZyEDKQ7hdK6mQdU8nMHgfK49z1FXd/ZLSHxbltyo/gxqoduBb4nLv/yMw+ROTsVddN\nZX2xEtSaBxQT+TX3XcDDZnZeJk+dmKDeLwPXT21Fo0vmZ9jMvkJkSuHfprK2JGTFZ2m8zGwO8CPg\nb9z9RKbricfMbgLa3P05M1uf7tfLunB394kEXjIn8U67sWo3sx8QmWsD2MwU/Fo2lgS1fhr4cTTM\nf29mA0T2wjg8VfWNNFq9ZnYBsAT4o5lB5L3/g5mtdffWKSxxSKKfYTP7BHATcG0Wnms4Kz5L42Fm\nM4gE+7+5+48zXc8YLgduNrMbgVnAPDP7P+7+0XS8WFimZbYAt5vZTDNbAiwHfp/hmkY6AFwVvXwN\n8FoGa0nkp0RqxMxWAPlk6aZM7v6Suy909xp3ryESTpdkKtgTMbN64IvAze7emel64ngWWG5mS8ws\nn0ijwpYM1zQqi/yL/j1gr7v/90zXMxZ3/5K7V0Z/Tm8HfpmuYIcsPHIfi5m9H/ifQBnwCzN7wd1v\niJ6w+2FgD5FfdT/j7v2ZrDWOTwH/YmZ5QDewMcP1jOUB4AEzawROA5/IwiPMoPo2MBN4LPqbxk53\n/4+ZLekMd+8zs7uAHUAu8IC7785wWWO5HPgY8JKZvRC97cvR8z5Pa1qhKiISQmGZlhERkRgKdxGR\nEFK4i4iEkMJdRCSEFO4iIiGkcBcRCSGFu4hICCncRURC6P8DlBR1FekPSx0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x3dd5518>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Beach erosion via scarp failure\n",
    "%matplotlib inline\n",
    "\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np \n",
    "\n",
    "# internal friction angle (is this the same as angle of repose?)\n",
    "phiw = 25.*(np.pi/180.)\n",
    "tanphiw = np.tan(phiw)\n",
    "# beach slope\n",
    "tanalphas = 1./50.\n",
    "alphas = np.arctan(tanalphas)\n",
    "# failure angle\n",
    "alpha = 0.25*np.pi+0.5*phiw*np.pi/180.\n",
    "tanalpha = np.tan(alpha)\n",
    "# dune slope\n",
    "tanalphad = 1./20.\n",
    "alphad = np.arctan(tanalphad)\n",
    "print(\"wet angle of repose: \",tanphiw,(180./np.pi)*phiw)\n",
    "print(\"beach slope: \",tanalphas,(180./np.pi)*alphas)\n",
    "print(\"failure angle: \",tanalpha,(180./np.pi)*alpha)\n",
    "print(\"dune slope: \",tanalphad,(180./np.pi)*alphad)\n",
    "\n",
    "xi = .2\n",
    "zs = 1.\n",
    "zi = xi*tanalpha\n",
    "# wedge of inundated sediment below slump block\n",
    "Vi = 0.5*xi*zi\n",
    "# wedge of dune on top of slump block\n",
    "Vd = 0.5*xi*xi*tanalphad\n",
    "# volume of slump block\n",
    "Ve = zs*xi-Vi+Vd\n",
    "print(\"Eroded volume: \",Ve,Vi,Vd)\n",
    "\n",
    "# first, try calcs assuming slump lies against slip face\n",
    "# we know the angles of the wedge of slumped sand\n",
    "A = phiw-alphas\n",
    "C = alpha-alphas\n",
    "B = np.pi-(A+C)\n",
    "print(\"A,B,C\",180./np.pi*np.array((A,B,C)))\n",
    "# and we know the area (volume per unit alonghore) of the wedge, so\n",
    "# use ASA to calculate sides\n",
    "# no cotangent function, use cot = 1/tan\n",
    "c = np.sqrt(Ve*2.*(1./np.tan(A)+1./np.tan(B)))\n",
    "b = np.sqrt(Ve*2.*(1./np.tan(A)+1./np.tan(C)))\n",
    "print(\"b, c = \",b,c)\n",
    "xs = c*np.cos(alphas)\n",
    "xl = b*np.cos(phiw)-xs\n",
    "print(\"xs, xl\",xs,xl)\n",
    "\n",
    "if xl > xi:\n",
    "    print('xl>xi...do calcs for slump resting against scarp')\n",
    "    xl = xi\n",
    "    Vis = 0.5*xi*xi*(tanalpha-tanphiw)\n",
    "    xs = np.sqrt(2.*(Ve+Vis)/(tanphiw-tanalphas))-xi\n",
    "    # check volume\n",
    "    Ve2 = 0.5*(xs+xl)**2.*tanphiw-0.5*(xs+xl)**2.*tanalphas-Vis\n",
    "    print(\"Ve, Ve2, diff:\", Ve,Ve2,Ve-Ve2)\n",
    "    print(\"xs, xl\",xs,xl)\n",
    "    \n",
    "# plot some points\n",
    "zb = .2\n",
    "zs = 1.\n",
    "xb = -zb/tanalphas\n",
    "zdt = zb+zs+.2\n",
    "xdt = (zdt-(zb+zs))/tanalphad\n",
    "zxs = zb-xs*tanalphas\n",
    "zslump = (xs+xl)*tanphiw+zxs\n",
    "x = np.array((xb,0.,0.,xdt))\n",
    "y = np.array((0,zb,zb+zs,zdt))\n",
    "plt.plot(x,y)\n",
    "plt.plot((0.,xi,xi),(zb,zb+zi,zb+zs+tanalphad*xi),'--')\n",
    "plt.plot((-xs,xl),(zxs,zslump),'-r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (<ipython-input-2-8f7877454480>, line 7)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  File \u001b[0;32m\"<ipython-input-2-8f7877454480>\"\u001b[0;36m, line \u001b[0;32m7\u001b[0m\n\u001b[0;31m    PHI =\u001b[0m\n\u001b[0m          ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "# q= -K(PHI)*dPHI/dx - eqn 3\n",
    "rho = 1025\n",
    "g = 9.81\n",
    "por = 0.3\n",
    "rhosb = por*rhos+(1.-por)*rho\n",
    "# PHI = pressure head (m)\n",
    "PHI = \n",
    "# K = hydraulic conductivity (m/s)\n",
    "Pw = rho*g*PHI\n",
    "# theta = volumentric water content of the pore spaces (m3/m3)\n",
    "# thetar = residual water content (initial value for a dry dune)\n",
    "# thetas = saturated water content\n",
    "# relative water content THETA i\n",
    "THETA = (theta-thetar)/(thetas-thetar)\n",
    "\n",
    "# Lateral earth force\n",
    "delz0 - zc-zb\n",
    "FE = rhosb*g*delz0/(2.*tanalpha2 )-Pw*THETA*(1.-1./tanalpha2)*delz0 # Eqn 21\n",
    "\n",
    "# Weight force of potential slump block\n",
    "FW = (rhosb*g+rho*g*theta)*(delzd-0.5*delxI*tanalpha)*delxI # Egn 23\n",
    "\n",
    "# Destabilizing force - eqn. 22\n",
    "FD = FW*np.sin(alpha)+FE*np.cos(alpha)\n",
    "FT = 2.*(-Pw)*THETA*tanphiw*np.tan(0.25*np.pi+0.5*phiw)*delzl\n",
    "# Stabilizing force - eqn. 25\n",
    "FS = (FW*np.cos(alpha)* + FE*np.sin(alpha) + (-Pw)*THETA*delxl/np.cos(alpha))*tanphiw+FT"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}