erodibility.ipynb

{
 "cells": [
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "MUD_TAUCR_MIN = 0.1\n",
    "MUD_TAUCR_MAX = 1.2\n",
    "MUD_TAUCR_SLOPE = 2.0\n",
    "MUD_TAUCR_OFF = -4.69e-1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "$$\\tau_{cb} = \\exp[(\\ln(M)-offset)/slope]$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ROMS/Nonlinear/Sediment/sed_bed_cohesive.F (COAWST v. 1085)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "!  Calculate representative critical shear stress profile\n",
    "!  Note that the values are for the TOP of the layer...we\n",
    "!  assume the bottom of the bottom layer has tcr = tcr_max\n",
    "          tcr(1) = tcr_min(ng)\n",
    "          DO k=2,Nbed\n",
    "             tcr(k) = tcr_min(ng)\n",
    "             IF (bmz(k-1).GT.eps) THEN\n",
    "                tcr(k) = exp((log(bmz(k-1))-                            &\n",
    "     &                    tcr_off(ng))/                                 &\n",
    "     &                    tcr_slp(ng))\n",
    "             ENDIF\n",
    "             tcr(k) = MIN( MAX( tcr(k), tcr_min(ng)), tcr_max(ng) )\n",
    "          ENDDO\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x8442a58>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAt8AAAJPCAYAAACkbLiHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XeUldXZxuH7mQJDG5AqvfcuRZo4oAiiYMsHktgTYzSW\naCxYEjAmwYK9xYgSNPaKgCI2FCmC9KaDUqT33pnZ3x/vEAlh+jlnn/K71poFzIwzdybL7c0++92P\nOecEAAAAIPySfAcAAAAAEgXlGwAAAIgQyjcAAAAQIZRvAAAAIEIo3wAAAECEUL4BAACACPFevs2s\nn5l9Z2aZZnZHLp/zhJktM7N5ZtYu0hkBAAHWbAAoHq/l28ySJD0lqa+klpKGmFmz4z7nbEkNnXON\nJV0j6R8RDwoAYM0GgBDwvfPdWdIy59wq59xhSa9LOu+4zzlP0kuS5Jz7RlJ5M6sW2ZgAALFmA0Cx\n+S7fNSWtPubPa3Lel9fnrD3B5wAAwo81GwCKyXf5BgAAABJGiufvv1ZSnWP+XCvnfcd/Tu18PkeS\nZGYupOkAIIKcc+Y7Qz5YswEgR1HXbN8737MkNTKzumZWQtLFkj447nM+kHSZJJlZF0k7nHMbc/uC\nzrlCv1X5Q38N+/f4Iv2zsfw2bNgw7xli6Y2fFz+vcL7FiKhYs/veN0Kn3n2H9//PIv3Gv1P8vPh5\nRc9bcXjd+XbOZZnZ9ZImKfiLwAvOuaVmdk3wYfdP59yHZtbfzH6QtFfSlT4zA0CiYs0GgOLzfexE\nzrmJkpoe977njvvz9RENBQA4IdZsACge38dO4FFGRobvCDGFn1fh8PMCQot/pwqHn1fh8POKHO87\n39GgzIFGKpWc7jtGxPEvWuHw8yocfl4Il7JWTVmujO8YEce/U4XDz6tw+HlFDjvfklr89LhalTvN\ndwwAQAF0SrlSe7+4QW+/LW3e7DsNABQOO98AgJhy+eVScrI0erT0619LdepIvXpJGRnS6adLlSr5\nTggAubPiXpcSTczMFeV/zznnSOedJ119tWTRfssugLhkZnLRf893SBV1zT7WkSPSnDnSF19IkydL\nU6dKDRoERbxXL6lnT+mkk0ISFwD+ozhrNuVb0quvSn/6U1C8/+//pEGDpHbtKOIAIofyHRqHD0vf\nfhsU8S++kKZPlxo3/nlnvGdPqXz5kH5LAAmI8p2jOAu5c8HuyVtvSW++GbykebSIt21LEQcQXpTv\n8Dh0SJo16+ed8W++kZo2Dcp4r15Sjx5SeuI9bw+gmCjfOYq6kC/bukxVy1RV+bRgO8Q5afbsn4t4\naurPRbxNG4o4gNCjfBfchj0bdCT7iGql1yr0P3vwYFDAj+6Mz5oltWz58854jx5S2bKF/rIAEgzl\nO0eRz3y/eo6u63idzmlyzv98zLngJcyjRbxkyZ+LeOvWFHEAoUH5Lrj7v75fOw7s0P1n3l/sDAcO\nSDNm/LwzPnt2sLYf3Rnv1k0qk3i3GgLIR3HWbK4azIeZ1KmT9OCD0ooV0ssvBzsnAwZIzZsHZ8UX\nLgxKOgAgMrJddki+TlpasON9773Sl19KmzZJf/1rcPTw3nulatWC3fBHH5W2bg3JtwSQ4CjfhWAm\nde4sPfSQtHKl9NJL0v790rnnBkX8z3+WFi2iiANAODU4qYEyt2aG5WuXLi2dcYZ0333S119LGzdK\n99wT7Ig3bCgNGSJ9/rmUHZruDyABUb6L6GgRHzkyKOJjxkh790r9+0stWkjDhkmLF/tOCQDxp0ut\nLpqxZoYicWyyTBmpXz/p3/+Wli8PjqH84Q9SkybSiBHS+vVhjwAgzlC+Q8BMOvVU6eGHgyI+erS0\ne3ewYLdsKQ0fLi1Z4jslAMSH2um1lWRJWrljZUS/b8WK0g03SPPnS6+8Iv34Y7DZcsEF0ocfSllZ\nEY0DIEZRviU1OqmR0kuG5q6ppCSpSxfpkUekVaukUaOknTulvn2DIn7vvdLSpSH5VgCQkMxMl7e9\nXNv2b/P0/YMNl1GjpJ9+ks4+O9hkqVcveNXzp5+8xAIQI7jtJEKys4Mn6t96S3rjDaljx+ChnjZt\nfCcDEA247ST2zZ8vPf+89NprwbHEq68OHs5PTfWdDECocdVgjlhZyA8ckJ57Ljgv2KtXsBvepInv\nVAB8onzHj337pLffDnbGMzOlK66Qfv3rYNImgPjAVYMxJi1Nuukm6YcfpFatpO7dgx2S1at9JwMA\nFFfp0tJll0lffRXcHX7kSLDO9+4tffwxN2IBiY7y7VHZstLddwc7I1WqSO3aBU/Rb9zoOxkAIBSa\nNQtuxVq9Otj9vukmqWfP4E5xAImJ8h0FTjpJ+vvfg6sJnQuenr/7bmn7dt/JACD6xcLRlZIlpV/9\nKpgFcfXV0lVXSWeeKU2f7jsZgEijfOeYuXamtyfnjzr5ZOnxx6W5c4Pd7yZNglK+Z4/XWAAQteZv\nmK9+r/QL2cTLcEtJCY6kfPedNHhw8HbOOcEQHwCJgfKd476v7tNXq77yHUOSVKdO8KDO1KnB6PrG\njaUnngjG2gMAftaqaivtOLBDY+aN8R2lUFJTgx3wZcuCqwoHDpQuvDBY8wHEN8p3jhaVW2jxpuga\nSdmkSXBl1UcfSZ98Evz5hReCh3cAAFJyUrKe6f+Mhn42VFv3bfUdp9BKlpSuvz54AL9Hj+AoysUX\nBzvjAOIT5TtHy6ottXhzdJXvo9q1k8aNk15/PZiq1qJF8Pvs2HiVFQDCqkONDhrUYpDu+uwu31GK\nrFQp6ZZbgqmZbdtKp50mXX55MNIeQHyhfOdoWaWl5qyfE9UP7nTtKn3+ufTss9Kjj0odOkhz5vhO\nBQD+3df7Po3LHKeZa2f6jlIsZctKd94Z7ITXrx8M67n2Wmlr7G3qA8gF5TtH++rtZWb6YuUXvqPk\n64wzgmmZt9wi9esn/eUv0uHDvlMBgD8V0ipozPljVL5ked9RQqJ8+WBk/fffB+fDW7QIjh3yiicQ\n+5hweYzxmeNVrkQ5nV7v9BCmCq81a6Tf/CbYFRkzJligAcQeJlwiL3PnBjvgSUnSM88ExxEB+MN4\n+RyJupA7Jz3/fHA3+NChwaCe5GTfqQAUBuUb+cnODna/77lHGjIkeNUzPd13KiAxMV4+wZlJv/2t\n9M030tixUkZG8NAOACB+JCUF1xMuXhzMf2jePHj4nr+/ALGF8h1HGjSQJk8O7oo99dTgwUwWZQCI\nL5UrB7Mg3npLGjFC6tMnOBsOIDZQvuNMUpJ0883SlCnS6NFS377S6tW+UwFAZK3csVJf//S17xhh\n1a1bMBnz3HODO8LvuUfat893KgD5oXznYvfB3dpzKHbnujdvLk2bJvXsKZ1ySvAwJrvgABLFTzt/\n0mXvXabDWfF9FVRKSvCcz/z5wXHDli2lSZN8pwKQF8p3Lm6ddKue+OYJ3zGKJSUl2An55BPp4Yel\nCy6QNm70nQoAwq9n3Z6qV6GeXl/0uu8oEVGjRjAR+bnnghuwrr+eXXAgWlG+c3F95+v15MwndeDI\nAd9Riq1dO2nWrOAawrZtpfff950IAMJvaI+hemDqA8p2iXM59llnSQsWSDt3Su3bSzNje+YQEJco\n37loXa21WldtrfeWvuc7SkiULCn9/e/Se+9JN94o/fnPDGsAEN/6NOijEsklNCFzgu8oEVWhgvTy\ny9Lf/iYNHCgNG8YgNiCaUL7zcPUpV+v5Oc/7jhFSXbsGu+BffCGdf760a5fvRAAQHmamoT2GauT0\nkb6jePGLXwTDeWbNCtb+pUt9JwIgUb7zdF6z87R482It27rMd5SQqlZN+uwzqVat4EpCrqgCEK8u\nan6RXrnwFd8xvKleXZowIbgfvGdP6fHHedUT8I0Jl/l4ZcEr6lCjg5pVbhbSrxstjk7G/Ne/pP79\nfacBEhcTLhFuP/wgXXaZVLp0cBVt7dq+EwGxi/HyOVjIi2bqVGnQoODp+KFDg4mZACKL8o1IOHJE\nevBB6bHHpCeflAYP9p0IiE2U7xws5EW3dm0wGbNu3WBHpEwZ34mAxEL5RiTNmRNsuvTrF1xFW7Kk\n70RAbCnOms2Zb0iSataUvvwyKN1du0orVvhOBAAIl1NOCaZjrlsnnXaatHKl70RA4qB84z/S0qQX\nXwwGNHTtGjyUCQDxYvGmxXFzfWwolC8vvfOOdPHFwcP3ExLrRkbAG8p3IRw8ctB3hLAzC+4Bf+01\n6Ve/Cs4F8qowgHhwKOuQbpx4o45kH/EdJWqYSbfcIr37rvS73wUP4B/hxwOEFeW7gFbvXK3mTzdP\nmElpvXpJM2YE579vuomrqQDEvvbV26tehXp6/zvG/B6ve/fgGMrMmcGUzA0bfCcC4hflu4Bql6+t\nciXLacqqKb6jREy9esE58PnzpUsukQ4d8p0IAIrnxs436olvnvAdIypVrSpNnBicAe/YUfrqK9+J\ngPhE+S6EX7b6pV5ZmFjDGipUCBbjvXuDMcV79/pOBABFd36z87VixwrNXT/Xd5SolJws3Xuv9MIL\nwW0oI0dy9BAINcp3IQxpPUTvLH0nIc5+H6tUqeChnBo1pDPOkLZu9Z0IAIomNTlV13S4RqPmjPId\nJar17RuMpX/tNemKK6QDB3wnAuIH5bsQ6pSvo1ZVW+njHz/2HSXiUlKCnZDTTw9ekly92nciACia\nGzrfoBFnjvAdI+rVri1NmSLt3y/17i1t3Og7ERAfKN+FdPUpV2vjnsRcgcykBx6QrrpK6tFD+u47\n34kAoPDKp5VXesl03zFiQunS0uuvBw9hdu4szZvnOxEQ+5hwiSIZM0a64w7pgw+CBRlA8TDhEtHu\nzTel3/9e+uc/pQsu8J0G8Ivx8jlYyCNr3LhgF/zVV6U+fXynAWIb5Rux4Ntvg+J97bXSnXcGr4gC\niYjynYOFPPK+/lq66CLpiSekwYN9pwFiF+UbsWLdOum886QmTaRRo4KH8oFEU5w1mzPfKJYePaRP\nPpH++MfggUwAiCXvLHlHew9xh2ph1KgRzIDIypIyMngQEygsyjeKrU0b6YsvpGHDpH/9y3caACi4\n0fNG652l7/iOEXNKlw6uIezXT+rWTVq2zHciIHZQvoto+fbluuuzu3zHiBqNG0uffSbdfbf00ku+\n0wBAwVze9nKNmT/Gd4yYZBYM5Bk6VOrZMxhNDyB/lO8iqlqmqp6e9bS279/uO0rUaNpU+vTT4CGc\nVxJrECiAGDWg6QDN3zBfq3as8h0lZl19dXADyjnnSB9+6DsNEP0o30VUtkRZ9a7fW+Myx/mOElWa\nN5cmTZJuvTW4GxYAollaSpoGtxyslxe87DtKTBswILh69qqrpBdf9J0GiG6U72K4qPlFnBU8gZYt\ngwJ+883SW2/5TgMAebui3RUaM3+MuHmleLp2DR7E/Otfpfvuk/hxAifGVYPFsOPADtV9rK7W3LxG\n5UqWi9j3jRXz50t9+0pPPx1cRwggd1w16I9zTh/98JH6Nuyr5KRk33Fi3vr1Uv/+0qmnBut/Mj9S\nxCGuGvSkQloFda/dXRN/mOg7SlRq21aaODGYiPbee77TAMCJmZn6N+5P8Q6R6tWDHfAffww2Xg4c\n8J0IiC7sfBfTxj0bVal0JaUkpUT0+8aSOXOks8+Wnn9eGjjQdxogOrHzjXhz6JB02WXS5s3S2LFS\n2bK+EwGhw4TLHCzk0evbb4OXIceMCYo4gP9G+UY8ysqSfvtbaenS4CaUChV8JwJCg2MniHodOwZP\nwl9+uTRjhu80AIBISE4OXvXs1Enq3TvYBQcSHeUbEdOlS7Dzff75wS4IAESbnQd2at/hfb5jxJWk\nJOmxx4JXPU8/XVq3znciwC/KNyLq7LOlBx8MRhKvXu07DQD8t6vHXa3XFzGkINTMpL/9LTgD3rOn\ntHKl70SAP5TvEFm3e5027d3kO0ZMuOwy6YYbggK+bZvvNADws8EtB+uVhYzoDZehQ6U//CEo4N9/\n7zsN4AflO0Qenvawnp31rO8YMePWW4MHMM89V9rHK7wAosQ5Tc7R3PVztXbXWt9R4tb110v33iv1\n6iUtXOg7DRB5lO8QGdB0AKPmC+mBB6RGjaRBg6TDh32nAYBg3PzApgP19pK3fUeJa1deKT3yiHTW\nWRRwJB7Kd4h0r91dy7cvZ7ekEJKSpBdekLKzg6uouHEMQDQY3HKw3lzypu8Yce/ii6VHHw0K+KJF\nvtMAkUP5DpHU5FT1a9RP4zPH+44SU1JTpbfeCm4/ufNO32kAQDqjwRlqXbW1jmQf8R0l7l188c87\n4IsX+04DRAblO4QGNOHoSVGUKSNNmBBMQHv0Ud9pACS6Eskl9I9z/8Hk4ggZMkQaOVLq00dassR3\nGiD8WFlCqF+jfpq3YZ7vGDGpUiXp44+lbt2kevWkCy7wnQgAECm//GVw9PDMM6VPP5VatPCdCAgf\nxssjqnz7bXAX+MSJUocOvtMAkcN4eUD697+l22+XPvtMat7cdxogd4yXR9zo2FF67rlgCuZanl0F\ngIRyySXBTVhnnil9953vNEB4cOwEUefCC6XMTGngQOmrr4Iz4QCAxHDppVJWVvAQ5ldfBUcRgXjC\nzjei0h13SK1bB7sg2dm+0wBIVJe8e4lW71ztO0bCueKK4PjJmWdK69b5TgOEFuUbUcksOH6ydat0\n112+0wBIVKnJqXrvu/d8x0hI118v/frXwS0oW7b4TgOEDuU7DBZvWqx7Pr/Hd4yYV7Kk9O670ttv\nS6NH+04DIBFd2OxCvbv0Xd8xEtaddwZHEPv1k3bu9J0GCA3KdxhULVNVT858UoeyDvmOEvMqV5bG\njw+OoXz5pe80ABJNn4Z9NG/DPG3au8l3lIT1979Lp54qnXuutG+f7zRA8VG+w6BKmSpqWqmppq2e\n5jtKXGjWTHr1VWnQIGnZMt9pACSStJQ09WvUT+9/977vKAnLTHrySal+/WAGxMGDvhMBxUP5DpN+\njfpp4g8TfceIG2eeKd13nzRggLRrl+80ABLJRc0vYj33LClJevFFqWzZYCDPkSO+EwFFx5CdMJmy\naopu/vhmffvbb31HiSvXXitt2CC9806wGAPxgiE70etw1mFJwcOX8OvgwWATpl694KF8S6h/YxBN\nGLIThU6tdaqWbVumrfu2+o4SVx5/XNq4MTgDCACRkJqcSvGOEiVLBpsvc+ZIw4f7TgMUDTvfYbRi\n+wrVrVBXScbfcUJp3TqpUydp1KhgFD0QD9j5Bgpu0yape3fplluCV0SBSCvOmk35RkyaOjWYhDl1\nqtSoke80QPFRvoHCWb5cOu006YknpIsu8p0GiSYmy7eZnSTpDUl1Ja2UNMg59z+3eJrZSkk7JWVL\nOuyc65zH12QhTyDPPCP94x/S9OmMoEfsi4XyHep1mzUbxTV3rtS3r/Tmm1JGhu80SCSxeuZ7qKRP\nnXNNJX0u6c5cPi9bUoZzrn1exRuJ59prpQ4dgglo/PcbiIiEXrd3HNihqT9N9R0Dx2jfXnr99eAq\n2vnzfacBCsZn+T5P0pic34+RdH4un2fiwVCcgJn07LPSjz9KjzziOw2QEBJ63d6wZ4MGvz1Y2S7b\ndxQco3dv6amnpHPOkVau9J0GyJ/PxbGqc26jJDnnNkiqmsvnOUmfmNksM7s6YulCxDnHZLQwSksL\nnnwfOVL67DPfaYC4lxDrdm6aVW6mMiXKaO76ub6j4DiDBkm33y717y9t3+47DZC3lHB+cTP7RFK1\nY9+lYFG+5wSfntvBge7OufVmVkXBYr7UOfd1iKOGzYY9G9Tq2VbadOsmJScl+44Tl+rUCSZgDhki\nzZol1a7tOxEQu1i383Zu43M1PnO8OtTo4DsKjnPjjdKKFcHD+B9/LJUo4TsRcGJhLd/OuT65fczM\nNppZNefcRjM7WdIJt4edc+tzft1sZu9J6iwp10V8+DEXf2ZkZCjD8xMY1ctVV5XSVTR/43ydUv0U\nr1niWa9ewZVTgwdLX34ppXIlL6Lc5MmTNXnyZN8x/kek1+1oW7PzM6DpAN32yW0aljHMdxScwMiR\n0i9+If3mN9KYMQzhQeiEcs32edvJA5K2OeceMLM7JJ3knBt63OeUlpTknNtjZmUkTZJ0r3NuUi5f\nMyqfnL92/LVqVLGR/tjtj76jxLXs7GDyWcuW0oMP+k4DFE6M3HYS0nU7WtfsvBzOOqxqI6tp8XWL\nVb1cdd9xcAL79gU3n/TvzyAehE+s3nbygKQ+Zva9pDMk3S9JZlbdzMbnfE41SV+b2VxJMySNy614\nR7Ne9Xvpi5Vf+I4R95KSpJdekt54Qxo/Pv/PB1BoCbNu5yY1OVUPn/WwjmQf8R0FuShdWho3Lvjv\nwZgx+X8+EGkM2YmATXs3qfGTjbX19q1KSQrrSR9ImjYtOPM3c2ZwHhyIBbGw8x1q0bpmIz4sXRrs\ngL/6qnTGGb7TIN7E6s53wqhapqr6NeqnjXs2+o6SELp1k/74x+D89+HDvtMAAHxo3jx4JXTIEGnx\nYt9pgJ+x8424lJ0tDRwYLL4PPeQ7DZA/dr6B8Hj55eDs9zffSJUr+06DeBGT4+XDgYUcx9q6VTrl\nlGD4woABvtMAeaN8A+Fz553S1KnSp59yBSFCg/Kdg4Ucx5s2TbrgguD8d926vtMAuaN8A+GTnS1d\ndJFUsaI0ahRXEKL4OPMN5KJbN+m22zj/DSD0Br89WJlbM33HQAEkJQXHT2bPlh591HcaJDrKN+Le\nLbdIlSpJw5iJASCE0kuk68NlH/qOgQIqW1b64INgEM+ECb7TIJFRviNo095NGj13tO8YCScpSRo9\nWvrXv6QoHCgIIEb1b9yf8h1j6tSR3nlHuvJKadEi32mQqCjfEZRsybpp4k0MZ/CgalXphRekyy6T\ntm/3nQZAPDijwRmavma69hza4zsKCqFr1+DoycCB0ubNvtMgEVG+I6hS6UqqW6Gu5q6f6ztKQjr7\n7ODhy2uukXjGC0BxpZdMV6canTR55WTfUVBIv/pV8CzQoEE8D4TIo3xHWM86PfXVqq98x0hYDzwQ\nTD371798JwEQD/o16qcvV37pOwaK4K9/lUqVkm691XcSJBquGoywNxe/qVcWvqKxF4/1HSVhLVok\n9eoVXEPYuLHvNECAqwZj04EjB1QiuYSSjL2sWLRjh9S5s3TXXdIVV/hOg1jCVYMx5LQ6p2nKqinK\ndtm+oySsVq2kP/85eNmRlxsBFEdaShrFO4ZVqCC9/35wJe3Mmb7TIFGwYkRY9XLV9WjfR3U4i9bn\n0/XXB2OGhw/3nQQA4FOLFtLzzwdDeDZs8J0GiYBjJ0hYGzdK7dpJr78unX667zRIdBw7AfwaPjwY\nP//554ygR/4YL5+DhRyF9eGH0rXXSgsWSOXL+06DREb5BvzKzg5uxKpeXfrHP3ynQbSjfOdgIUdR\nXHNNcPb7xRd9J0Eio3zHtlU7Vik5KVm10mv5joJi2LUreADzjjuCQTxAbnjgEiiGkSODyZfjxvlO\nAiBWPT/neT0z6xnfMVBM6enSu+9Kt98uzWUkB8KE8o2EV65ccO/3NddIW7b4TgMgFp3V8CxN+nGS\n7xgIgRYtpKefDh7A3LbNdxrEI8q3Jx8t+0h3f3a37xjI0bOnNGSIdN11TL8EUHhdanXRsm3LtGUf\nf4OPB4MGSeedJ116aXAWHAglyrcn1cpW03vfvec7Bo7xt78FA3hef913EgCxpkRyCfWs21OfLf/M\ndxSEyIMPBmfA//pX30kQbyjfnrSp1kZrdq3R1n1bfUdBjrQ06aWXpJtuktat850GQKw5q8FZ+mT5\nJ75jIERSU6U335See06aONF3GsQTyrcnKUkp6lKri6atnuY7Co7RsWNw9eBvfsPxEwCF079xfzWt\n1NR3DIRQ9erBq6GXXy6tXOk7DeIF5dujHnV66OufvvYdA8e5555gAM+oUb6TAIglDSs21G3db/Md\nAyF22mnB1YODBkmHDvlOg3hA+faoR50emrp6qu8YOE5qanD85K672OkAAEg33yzVqCHdxt+tEAIM\n2fHoUNYhHThyQOkl031HwQncf38wZvjjjyVLqNEn8IEhO0B0275d6tBBeuih4BpCJDYmXOZgIUco\nHTkidekSXD941VW+0yDeUb6B6DdrlnTOOdL06VLDhr7TwCfKdw4WcoTa/PlSnz7SvHnBS45AuFC+\ngdjw5JPS6NHStGnBLVlITJTvHCzkCIdhw4Ly/f77HD9B+FC+48eEzAnatn+bLm17qe8oCAPngocv\nq1SRnnnGdxr4Upw1mwcugXzcfbf044/SG2/4TgIgFmS7bI2ZP8Z3DISJWXAb1qRJDGVD0bDzHQX2\nHNqjJEtS6dTSvqMgFzNnSgMHSgsXBrsdQKix8x0/dh3cpZqP1NTm2zYrLYVzCfFqzhypb19pxgzO\nfycidr5j3OXvX66x3431HQN56NxZuvRS6cYbfScBEO3SS6arVdVWmr56uu8oCKNTTpH+9Cfp4ou5\n/xuFQ/mOAl1rdWXSZQz4y1+k2bOlsfw9CUA+etfrrc9XfO47BsLshhuCKZh33eU7CWIJ5TsKdK3V\nVdPWUL6jXalSwTm/664L7nsFgNz0rt9bn6+kfMc7M+nFF4Nngj76yHcaxArOfEeBA0cOqOIDFbXl\n9i2c+44B110nZWVJzz3nOwniCWe+48uBIwe0csdKNavczHcURMCXXwbHT+bMCXbCEf848x3j0lLS\n1KpqK81ZP8d3FBTAiBHShAnS11/7TgIgWqWlpFG8E8jpp0vXXCNdckmwOQPkhfIdJQY0GaCt+7b6\njoECKF9eevxx6be/lQ4e9J0GABAN7rknmIz84IO+kyDacewEKALnpPPOkzp1Cp52B4qLYydA7Js9\nO7gZa8kS30kQbhw7ASLMTHrqqWAHPDPTdxoAQDQoVy7Y/QbyQvkGiqhOneBlxt/9LtgJB4DjOee0\n7/A+3zEQIcnJlG/kj/INFMMNN0i7dkljmCQN4ARGThupP33O2bREkZLCA5fIH+UbKIbkZOn556U7\n7pA2b/adBkC06Va7myavmuw7BiKEnW8UBOU7iuw+uFvvLn3XdwwUUvv2wQM2t9ziOwmAaNOpZidl\nbs3U9v31lgHqAAAgAElEQVRM5koE7HyjICjfUcTJ6bL3LtPhrMO+o6CQ7r1XmjJF+pyBdgCOUSK5\nhLrW6qopP03xHQURkJxM+Ub+KN9RJL1kuupWqKvFmxf7joJCKlMmuPnk97+XDh3ynQZANMmol6HJ\nKyf7joEISEqifCN/lO8o07lmZ81cO9N3DBTBwIFSo0bSI4/4TgIgmvSq10s7DuzwHQMRYAl1Uz+K\niiE7UebZWc9q9vrZGjVwlO8oKIIVK4LBO7NnS3Xr+k6DWMKQHSD2bdsWbMJs2+Y7CcKNITtxpHPN\nzvpm7Te+Y6CI6teX/vCH4A0AAOB4lO8o07paa13U/CKxGxS7brtNWrxYmjDBdxIAABBtOHYChMGk\nScHky8WLpVKlfKdBLODYCRD7OHaSODh2AkSZs86SOnaURozwnQQAAEQTyjcQJo88Ij3zjLRsme8k\nAKLB9NXT9f2W733HAOAZ5RsIk1q1pDvvlK6/XuKVdQDjM8fr3wv+7TsGAM8o30AY3XijtHq19MEH\nvpMA8K1n3Z766qevfMcA4BnlO0o9O+tZTV893XcMFFNqajD58uabpQMHfKcB4FPX2l01e91sHTxy\n0HcUAB5RvqPUTzt/0qQfJ/mOgRDo00dq1056+GHfSQD4lF4yXc0qN9OsdbN8RwHgEeU7SnWs0VHf\nrv/WdwyEyCOPSI8+GhxBAZC4TqtzmqasmuI7BgCPKN9RqlPNTvp23bcM24kT9epJv/99MIAHQOL6\nZetfqn319r5jAPAoxXcAnFjt9NrKys7S2t1rVSu9lu84CIE77pCaN5e+/FI6/XTfaQD40KlmJ98R\nAHjGzneUMrP/7H4jPpQuHZz7vuEG6cgR32kAAIAPlO8odv8Z96trra6+YyCELrpIqlxZeu4530kA\nAIAPFk9nis3MxdP/HsSnRYuk3r2lJUuCIg5IwatdzjnznSOSWLMRb7Ztkxo1Cn5FfCvOms3ONxBh\nrVpJF18sDR/uOwkAAIg0yjfgwbBh0ptvBrvfABLPteOv1aodq3zHAOAB5RvwoFIl6e67pT/+0XcS\nAD5s2b9FU37ivm8gEVG+AU+uu05asUL66CPfSQBEWo/aPfT1T1/7jgHAA8p3lJu+eroufvti3zEQ\nBqmpwdWDt9wiHT7sOw2ASOpep7umrp7qOwYADyjfUa5uhbr6dPmnTLqMU/37S7VrS//4h+8kACKp\nbbW2Wrljpbbv3+47CoAIo3xHuRrlaig1OVWrdvJgTjwykx55RLrvPq6mAhJJanKqOtXopOlrpvuO\nAiDCKN8xoEP1Dpq9brbvGAiTVq2kX/xC+stffCcBEEkvnveiMupl+I4BIMIo3zGgQ/UOmr2e8h3P\n7r1XeuUV6bvvfCcBECn1KtRT6dTSvmMAiDDKdwzoUKOD5m2Y5zsGwqhKFWnoUOn2230nAQAA4cR4\n+RhwKOuQJKlEcgnPSRBOBw9KzZpJo0dLGRm+0yDSGC8PxD7GyycOxsvHuRLJJSjeCaBkSWnECOnW\nW6XsbN9pAABAOFC+gSgyeLCUlCS9/rrvJAAi5eCRg74jAIggyjcQRcykkSOlu+6SDhzwnQZAuK3a\nsUpNnmriOwaACKJ8A1GmZ0+pXTvpiSd8JwEQbnXK19H+w/u1eudq31EARAjlO4Zs2rvJdwREyAMP\nSA89JG3Z4jsJgHAyM3Wt3ZVhO0ACoXzHCOecmjzZRJv3bvYdBRHQtKk0aFAw+RJAfOtWq5umrZ7m\nOwaACKF8xwgzU/vq7TV3w1zfURAhw4YFg3eWLfOdBEA4datN+QYSCeU7hrQ/ub3mrJ/jOwYipGpV\n6dprpaef9p0EQDh1rNFRew7tUVZ2lu8oACKA8h1D2p/cnkmXCaZOHWnPHt8pAIRTqdRSWvL7JUpO\nSvYdBUAEUL5jCMdOEk9yspTFZhgAAHGD8h1DmlVupvSS6Tqcddh3FERIUhLTLgEAiCcpvgOg4FKS\nUjTr6lm+YyCC2PkGACC+eNv5NrNfmNkiM8sys1Py+Lx+ZvadmWWa2R2RzAj4xs43ognrNgAUn89j\nJwslXSDpy9w+wcySJD0lqa+klpKGmFmzyMQD/GPnG1GGdTuMMrdmKnNrpu8YAMLMW/l2zn3vnFsm\nyfL4tM6SljnnVjnnDkt6XdJ5EQkIRAF2vhFNWLfD672l7+nZWc/6jgEgzKL9gcuaklYf8+c1Oe8D\nEgI734hBrNtF1KVWF81YO8N3DABhFtYHLs3sE0nVjn2XJCfpbufcuHB+73j2zZpvVK9CPVUrWy3/\nT0ZMS0qifCOyWLf96VijoxZsXKCDRw6qZEpJ33EAhElYy7dzrk8xv8RaSXWO+XOtnPflavjw4f/5\nfUZGhjIyMooZIfo8OuNRndvkXF3S5hLfURBmSdH+2hSKbPLkyZo8ebLvGP8j0ut2IqzZBVWmRBk1\nqdRE8zbM06m1TvUdB8AxQrlmm3MuJF+oyAHMvpB0q3Nu9gk+lizpe0lnSFovaaakIc65pbl8Lef7\nf08kjJgyQlv3b9XIs0b6joIwGztWevHF4FfENzOTcy6vs9RRI1TrdqKs2YVx7fhr1axyM93U5Sbf\nUVAE27ZJjRoFvyK+FWfN9nnV4PlmtlpSF0njzeyjnPdXN7PxkuScy5J0vaRJkhZLej234p1I2p3c\njjHzACKOdTv8Lmx+oWqmc0QeiGfed75DKVF2UdbtXqe2/2irTbdukllMbJShiNj5ThyxtPMdKomy\nZiNxsPOdOGJy5xtFV71sdZlM63av8x0FAAAAhUD5jkFmpptOvUkHsw76jgIAAIBCCOttJwifu3ve\n7TsCAAAAComdbwAAACBCKN8AAESRI9lHdMvHtyjbZfuOAiAMKN8AAESRlKQUvf/d+8rcmuk7CoAw\noHwDABBlTq11qr5Z843vGADCgPIdwz758RNNWTXFdwwAQIh1rtFZM9fO9B0DQBhQvmPYgo0L9NaS\nt3zHAACEWOeanTVzHeUbiEeU7xjW9uS2WrBxge8YAIAQa1+9vRZvWqwDRw74jgIgxCjfMaxNtTaa\nv3G+GM8MAPGldGppvfV/vLIJxCPKdwyrWqaq0lLStGbXGt9RAAAhdk6Tc5SWkuY7BoAQo3zHuKO7\n3wAAAIh+jJePcbd1u02102v7jgEAAIACoHzHuDMbnOk7AgAAAAqIYycAAABAhFC+AQCIUs/OelZP\nfvOk7xgAQojyDQBAlKqQVkFfrPzCdwwAIUT5BgAgSnWq2Umz1s3yHQNACFG+48DEHybq6ZlP+44B\nAAixhic11N5De7VhzwbfUQCECOU7DmRlZ+mDzA98xwAAhJiZBbvfa9n9BuIF5TsOtK7WWgs2LvAd\nAwAQBh2rd9S36771HQNAiHDPdxyonV5b+w/v1+a9m1WlTBXfcQAAITS0x1DGzANxhJ3vOGBmal2t\ntRZuWug7CgAgxMqVLKfU5FTfMQCECOU7TrSu2lqLNi3yHQMAAAB54NhJnLi12628LAkAABDlKN9x\nosFJDXxHAAAAQD44dgIAQAzYvn+77wgAQoDyDQBAlNu8d7MaPNFAzjnfUQAUE+UbAIAoV6VMFZUr\nUU7Lty/3HQVAMVG+AQCIAR1qdGDYDhAHKN9x5MNlH+p343/nOwYAIAw6Vu+o2etn+44BoJgo33Hk\n5LIna+rqqb5jAADCgJ1vID5QvuNI88rN9cO2H3Q467DvKACAEOtQvYMOZ7O+A7GO8h1HSqWWUu30\n2lq2bZnvKACAEKtSpoqmXDnFdwwAxUT5jjOtqrZizDwAAECUynPCpZmlSTpX0mmSakjaL2mRpAnO\nucXhj4fCalW1lRZvWiy19J0EAAAAx8u1fJvZvQqK92RJ30jaJClNUhNJ9+cU8z865xZEICcK6I7u\nd6hkSknfMQAAAHACee18z3TODcvlY4+YWVVJdcKQCcVQpkQZ3xEAAACQi1zPfDvnJuT1DzrnNjnn\nuPMIAIAImrZ6mrbs2+I7BoAiyrV8m1mymV1jZveZWffjPnZP+KMBAIDjPTTtIX22/DPfMQAUUV63\nnTwn6XRJWyU9YWaPHPOxC8OaCgAAnFCH6h2YdAnEsLzKd2fn3C+dc49JOlVSWTN718xKSrLIxENR\nOOe099Be3zEAAGHQoXoHzd0w13cMAEWUV/kucfQ3zrkjzrnfSpon6XNJZcMdDEU39vuxGvLOEN8x\nAABhcEr1UzR73Ww553xHAVAEeZXvb82s37HvcM79RdJoSfXCGQrF07xycy3ezDXsABCPqpWtplKp\npbRq5yrfUQAUQV63nVzinJt4gvePcs6lhjcWiqNhxYZav3u99h3e5zsKACAMbjr1Jh3KOuQ7BoAi\nyHe8vJklRyIIQiclKUWNKzXW0s1LfUcBAITB7d1vV5NKTXzHAFAEeZZvMysnaWyEsiCEWlZpydET\nAACAKJPXPd/VJX0q6Z+Ri4NQOaX6KQxhAAAAiDJ5jZefIuk259wHkQqD0Lm1262+IwAAAOA4eR07\n2S6pZqSCAAAAAPEur/KdIelsM/t9hLIAAIACmrZ6mj74nhengViT11WDeyUNlNQ+cnEAAEBBrNm1\nRqPnjfYdA0Ah5XnbiXMuyzn3m0iFAQAABdP+5Paau54x80Csyfee76PMLN3MKh59C2cohMb2/dv1\nw7YffMcAAIRBw4oNtW3/Nm3dt9V3FACFUJAhO9eY2QZJCyTNznn7NtzBUHxfrPxCN398s+8YAIAw\nSLIktTu5neZuYPcbiCUF2fm+VVIr51w951z9nLcG4Q6G4mtZpaWWbF7iOwYAIEw4egLEnrzu+T7q\nR0n7wh0EodewYkOt271O+w7vU+nU0r7jAABC7Mr2VyorO8t3DACFUJDyfaekaWb2jaSDR9/pnLsx\nbKkQEilJKWpUsZG+3/K92lfn0hoAiDftTm7nOwKAQipI+X5O0ueSFkrKDm8chFqLKi20ZPMSyjcA\nAEAUKEj5TnXO3RL2JAiLPg36KCWpIP83A4h1ZpYm6VxJp0mqIWm/pEWSJjjnFvvMBgAIFKSVfWRm\nv5U0Tv997GRb2FIhZH5zCte0A4nAzO5VULwnS/pG0iZJaZKaSLo/p5j/0Tm3wFtIAECByveQnF/v\nPOZ9ThI3ngBA9JjpnBuWy8ceMbOqkupEMhAA4H/letWgmdWQpGOuF6zPVYMAEJ2ccxPy+fgm5xwz\nGuLQ8u3LddXYq3zHAFBAed3zPcrMZpjZ/WaWYWYcHAaAKGdmHc3sPTObY2YLzGyhmXHUJI5VKV1F\nbyx+Q0eyj/iOAqAAci3Uzrn+OWcEMyRdIGmkmf0kaaKkic65nyITEQBQCK9Iuk3cUJUwypUsp5rl\naipza6ZaVGnhOw6AfOS5m+2cO6Ccsi1JZlZf0tmSnjKzk51zncMfEcWVuTVT2/ZvU5daXXxHARB+\nm51zH/gOgchqd3I7zdswj/INxICCjJf/D+fcCufcM865gZJ6hCkTQmzm2pl6bMZjvmMAiIxhZjbK\nzIaY2YVH33yHQni1O7kdY+aBGJFv+Taz3Wa267i31ZLeMDMevIwBLaq00NItS33HABAZV0pqJ6mf\npAE5b+d6TYSwa39ye83bOM93DAAFUJCHKB+TtEbSq5JM0sWSGkqaI+lFBWfCEcWaVmqqzK2ZysrO\nUnJSsu84AMKrk3Ouqe8QiKyedXuqWeVmvmMAKICCHDsZ6Jx7zjm32zm3yzn3T0l9nXNvSDopzPkQ\nAmVKlFG1MtW0fPty31EAhN80M+Pgb4IpU6KM6p9U33cMAAVQkPK9z8wGmVlSztsgSQdyPubCmA0h\n1LxKc3235TvfMQCEXxdJ88zse64aBIDoU5BjJ7+S9LikZxSU7RmSLjGzUpKuD2M2hNClbS5V5dKV\nfccAEH79fAcAAOSuIOV7h3NuwLHvMLP6zrn9kr4OTyyE2i9b/9J3BAAR4JxbZWYnSaqt/17jV3mK\nBAA4RkHK9zgzO9s5t0uSzKy5pLcktQprMgBAoZnZfZKukPSjfj4a6CT19pUJkeOck5n5jgEgDwUp\n339XUMDPkdRU0ksKjqIAAKLPIEkNnXOHfAdBZGW7bNV8pKaW3bBMZUuU9R0HQC7yLd/OuQlmlipp\nkqRyki5wzmWGPRkAoCgWSaogaZPvIIisJEtSrfRaWrBxgbrV7uY7DoBc5Fq+zexJ/fdtJuUVvIx5\nvZnJOXdjuMMBAApthKS5ZrZI0sGj78yZTIw417ZaW83fMJ/yDUSxvHa+vz3uz7PDGQTh9+nyT5WS\nlKKMehm+owAInzGSHpC0UFK25yyIsLbV2mr+xvm+YwDIQ67l2zk3JpJBEH4LNy7U8u3LKd9AfNvn\nnHvCdwj40fbktnp10au+YwDIQ65DdsxsnJkNyDnvffzHGpjZX8zsqvDGQyg1q9xM321l0A4Q56aY\n2Qgz62pmpxx98x0KkdGmWhst375czjEDD4hWeR07uVrSLZIeM7NtkjZLSpNUX9IPkp5yzo0Nf0SE\nClMugYTQPufXLse8j6sGE0SFtApad8s6rhsEolhex042SLpd0u1mVk9SdUn7JWU65/ZFJB1Cqk75\nOtq6b6t2H9ytciXL+Y4DIAycc718Z4BfyUnJviMAyEOux06O5Zxb6Zyb7pybR/GOXUmWpCaVmuj7\nrd/7jgIgxMzsEjPL6yhhQzPrEclMAID/VZAhO2FhZr+QNFxSc0mdnHNzcvm8lZJ2Knhq/7BzrnOk\nMsaj4RnDVbVMVd8xAIReJQVXDM5WcDvV0aOCjSSdLmmLpKHF+Qas2wBQfN7Kt4JrsC6Q9Fw+n5ct\nKcM5tz38keLf+c3O9x0BQBg45x43s6cUnO3uLqmNgqOCSyVd6pz7KQTfhnUbAIrJW/l2zn0vSZb/\nUyGmAh6PAYBE5pzLkvRJzls4vj7rdoxYvXO1qpSporSUNN9RABwn38XRzLqb2Sdmlmlmy81shZkt\nj0S4HE7SJ2Y2y8yujuD3BQAUDeu2Z4PfHqyZa2f6jgHgBAqy8/2CpJsVnCHMKswXN7NPJFU79l0K\nFuW7nXPjCvhlujvn1ptZFQWL+VLn3NeFyQEAKBjW7fjQplobLdi4QD3r9vQdBcBxClK+dzrnPirK\nF3fO9SnKP3fc11if8+tmM3tPUmdJuS7iw4cP/8/vMzIylJGRUdwIABBykydP1uTJk33H+B+RXrdZ\ns8OjTbU2mrt+ru8YQNwI5ZptuU3BOmYi2iBJyZLelXTw6Mdze8q90AHMvpB0q3Nu9gk+VlpSknNu\nj5mVkTRJ0r3OuUm5fC3HVK/8PTbjMXWu2VndanfzHQX5GDtWevHF4FfENzOTc67Yk1HMrKSkiyTV\n0zEbLM65vxT3ax/zPUKybrNmh8/Un6bq5o9v1syrOXoSSdu2SY0aBb8ivhVnzc5r5/vh4/7c8Zjf\nF3tampmdL+lJSZUljTezec65s82suqTnnXPnKnjp8z0zczlZX8mteKPgVu5YqazsLMo3EJ/GKrjm\nb7aO2TAJBdbt2NGqaist2bxEWdlZDN0BokxeEy57SZKZNXDO/dcDlmbWoLjf2Dn3vqT3T/D+9ZLO\nzfn9Ckntivu98N+aVmqqOetD8sIFgOhTyznXLxxfmHU7dpRPK69e9Xtp6/6tzHYAokxBroJ6+wTv\neyvUQRA5zSo3Y8olEL+mmVlr3yHg37gh4yjeQBTKdefbzJpJaimpvJldeMyH0hVMTUOMalq5qb7b\n8p3vGABCyMwWKjgSmCLpypwrYQ8q57YS51wbn/kAAIG8znw3VfAyYgVJA455/25J3Nsaw6qXra4D\nRw5o+/7tOqnUSb7jAAiNc30HAADkL68z32MljTWzrs656RHMhDAzM30w5AMmnwFxxDm3SpLM7GXn\n3KXHfszMXpZ06Qn/QQBARBXkzPdGMxtnZpvNbJOZjQ3FA5fwK6NehkqllvIdA0DotTz2D2aWLKmD\npywAgOMUpHy/KulNSdUl1VDwsOVr4QwFACgcM7vTzHZLamNmu8xsd86fNym4fhAJaPrq6Vq9c7Xv\nGACOUZDyXdo597Jz7kjO27/FA5cAEFWccyOcc+UkPeScS3fOlct5q+Scu9N3Pvgxas4ofbjsQ98x\nAByjIOPlPzKzoZJeV/Ak/WBJH5pZRUlyzjHHCQCix105N1T1ULBmT8m5nxsJqHW11lqwcYHvGACO\nUZDyPSjn12uOe//FChZ2zn8DQPR4WlIj/Xw88Hdm1sc593uPmeBJm2pt9O7Sd33HAHCMfMu3c65+\nJIIg8i568yI91OchNTiJvz8BcaS3pObOOSdJZjZG0mK/keBL66rBzrdzTmbmOw4AFeDMt5mVNrN7\nzOyfOX9ubGbcJxsH9hzao6Wbl/qOASC0fpBU55g/1855HxJQlTJVVCq1lFbv4qFLIFoU5IHL0ZIO\nSeqW8+e1kv4atkSImKaVmjJmHog/5SQtNbPJZvaFpCWS0s3sAzP7wHM2eHBzl5t1JPuI7xgAchTk\nzHdD59xgMxsiSc65fcZrV3GhaaWmPIgDxJ8/+w6A6HJ799t9RwBwjIKU70NmVkrBw5Uys4aSDoY1\nFSKiaeWmemvJW75jAAgh59yXZlZXUmPn3Kc563eKc26372wAgIIdOxkmaaKk2mb2iqTPJPHX6DjA\nsRMg/pjZ1ZLelvRczrtqSeKqQQCIEgW57eQTM5sjqYskk3STc25L2JMh7Gqm19TXV37tOwaA0Pq9\npM6SvpEk59wyM6vqNxIA4Kg8y7eZpUg6W1KznHctlbQj3KEQGUmWpIYVG/qOASC0DjrnDh19NCdn\nHXd+IwEAjsr12ImZ1VRwN+wfJdWQVFPSbZIWm1mNyMQDABTSl2Z2l6RSZtZH0luSxnnOBM/GZ47X\n9NXTfccAoLx3vv8m6Vnn3GPHvtPMbpQ0QtLl4QwGACiSoZJ+LWmhgsnEH0oa5TURvJuzfo72H96v\nrrW7+o4CJLy8yncX59wVx7/TOfeEmfGUHgBEIedctpm9L+l959xm33kQHVpXba0X573oOwYA5X3b\nyf48PrYv1EEAAEVngeFmtkXS95K+N7PNZsa931Draq21cONC3zEAKO+d7/JmduEJ3m+S0sOUBxG2\n59AeNX2qqdbcvEbMTgJi2s2Sukvq5JxbIUlm1kDSs2Z2s3PuUa/p4FX9CvW1ed9m7Tq4S+kl+U84\n4FNe5ftLSQNy+dhXYcgCD8qWKKvDWYe1Yc8GVS9X3XccAEV3qaQ+x14F65xbbmaXSJokifKdwJKT\nktWiSgst3rSYc9+AZ7mWb+fclZEMAn+aVg6G7VC+gZiWeqIZDM65zWaW6iMQosu9GfeqRjkuKwN8\nK8iES8S5JhWbKHNrpu8YAIrnUBE/hgTRv3F/1a1Q13cMIOHlO+ES8a9JJco3EAfamtmuE7zfJKVF\nOgwA4MQo31DjSo01Y+0M3zEAFINzLtl3BgBA/vI9dmJm/2dm5XJ+f4+ZvWtmp4Q/GiJlYNOBemfQ\nO75jAAAAxL2CnPn+k3Nut5n1kHSmpBckPRveWIiklKQUJRnH/wEAAMKtII0rK+fXcyT90zk3QVKJ\n8EUCAADh8O8F/9bri173HQNIaAUp32vN7DlJgyV9aGYlC/jPAQCAKLLn0B59uvxT3zGAhFaQEj1I\n0seS+jrndkiqKOm2sKYCAAAh16pqKy3atMh3DCChFaR8V5c0wTm3zMwyJP2fpJlhTYWIy3bZ2nlg\np+8YAIAwalmlpRZvXqxsl+07CpCwClK+35GUZWaNJP1TUm1Jr4Y1FSLu4x8+1i/e+oXvGACAMDqp\n1EkqX7K8ftr5k+8oQMIqSPnOds4dkXShpCedc7cp2A1HHGHQDgAkBo6eAH4VZMjOYTMbIukySQNy\n3pcavkjwoW6Futq4Z6MOHDmgtBSG4QFAvHq076OqVraa7xhAwirIzveVkrpK+ptzboWZ1Zf0cnhj\nIdJSklJUt0Jd/bjtR99RAABh1LxKc1UsVdF3DCBh5Vu+nXNLnHM3Oudey/nzCufcA+GPhkjj6AkA\nAEB45XvsxMwaSxohqYWk/5xHcM41CGMueNCuWjttP7DddwwAAIC4VZAz36MlDZP0qKReCo6hMGQn\nDt3X+z7fEQAAAOJaQUp0KefcZ5LMObfKOTdcwah5AAAAAIVQkJ3vg2aWJGmZmV0vaa2ksuGNBQAA\nwuWpmU/pwJEDurXbrb6jAAmnIDvfN0kqLelGSR0kXSrp8nCGAgAA4VMhrYK+Xfet7xhAQirIbSez\nnHN7nHNrnHNXOucudM7NiEQ4AAAQegzaAfzJ9diJmX2Q1z/onBsY+jjwbf3u9UqyJAYwAEAca1a5\nmX7c/qMOZx1WajJz84BIyuvMd1dJqyW9JukbSRaRRPDqsRmPqXxaed112l2+owAAwiQtJU2102tr\n2bZlalGlhe84QELJ69jJyZLuktRK0uOS+kja4pz70jn3ZSTCIfIaVWykZduW+Y4BAAizllVbavGm\nxb5jAAkn1/LtnMtyzk10zl0uqYukHyRNzrnxBHGqcaXGWraV8g0A8W7M+WN0UYuLfMcAEk6eVw2a\nWUkFd3oPkVRP0hOS3gt/LPjSuGJjdr4BIAGkl0z3HQFISHk9cPmSgiMnH0q61znHY9EJoEa5Gtpz\naI92HdzFwgwAABBieZ35vkRSYwX3fE8zs105b7vNbFdk4iHSzEznNztf2/Zv8x0FAAAg7uS68+2c\nK8gAHsShVy58xXcEAACAuETBBgAgQWW7bO0/vN93DCChUL4BAEhQ9315n0Z8PcJ3DCChUL4BAEhQ\nTSs31ZLNS3zHABIK5RsAgATVskpLLd7MoB0gkijfOKFJP07SvsP7fMcAAIRRk0pNtGL7Ch3KOuQ7\nCpAwKN84ods+uU1LNy/1HQMAEEYlU0qqXoV6ytya6TsKkDAo3zihRhUb6YdtP/iOAQAIs041O2nt\nrrW+YwAJI8/x8khcjSs2pnwDQAJ4+YKXfUcAEgo73zihRhUbadm2Zb5jAAAAxBXKN06InW8AAAqn\nXLH+9JsAAA9wSURBVDnp4499p0C0o3zjhJpVbqbONTv7jgEAQMxITZU6dfKdAtGO8o0Tqla2mh7p\n+4jvGAAAAHGF8g0AAABECOUbAAAAiBDKNwAAABAhlG8AAAAgQijfyNXmvZv12sLXfMcAAACIG5Rv\n5Grf4X267ZPbfMcAAACIG5Rv5KpWei1t2bdF+w7v8x0FAAAgLlC+kavkpGTVq1BPK7av8B0FAAAg\nLlC+kadGFRsxZh4AACBEKN/IU8OTGurH7T/6jgEAABAXUnwHQHQ7r9l5ysrO8h0DAAAgLlC+kafe\n9Xv7jgAAABA3OHYCAAAARAjlGwAAAIgQyjcAAAAQIZRvAAAAIEIo38jXR8s+0mfLP/MdAwAAIOZR\nvpGvRZsWaXzmeN8xAAAAYh7lG/lqWJFBOwAAAKFA+Ua+mHIJAAAQGt7Kt5k9aGZLzWyemb1jZum5\nfF4/M/vOzDLN7I5I50Sw871i+wplu2zfUQB4xLoNAMXnc+d7kqSWzrl2kpZJuvP4TzCzJElPSeor\nqaWkIWbWLKIpobIlyiq9ZLrW717vOwoAv1i3AaCYvJVv5/6/vbsNsfQs7wD+v7KTuMbQzMYUI4lJ\nZmNVmhYTxLTFkkwVaRSD/SBUKYqKxbao/SL0JWK2pR/SgGBf8ENEkRbfoGCNL5FsKPNBMKmoqWJt\nsu4mUUMqsmZtY8Suyd0PM8bNZs7unDln7uc8s78fDHOec27OXHNzn2v/+8xzzt3ubO3JU6l3Jblk\nk2HXJDnUWnuwtXY8ySeSvLZXjfzC+69/f/Yu7R26DGBA+jbA7Bblmu+3Jrl9k/svTvLdE46/t3Ef\nnb3+116fZ5/77KHLABaHvg2wDUs7+eRVdTDJc068K0lLcmNr7TMbY25Mcry19rGdrAWA09O3AXbW\njobv1torT/V4Vb05yauTvHzCkIeSXHrC8SUb90104MCBJ2+vrq5mdXX19IUCdLa2tpa1tbWhy3ia\n3n1bzwbGYJ49u1prc3miqX9w1fVJ3pfk2tba0Qlj9iS5N8krkjyc5N+TvKG19q0J49tQvw/shE9/\nOvnwh9e/s7tVVVprNXQdpzLvvq1nA2M1S88e8prvf0hyXpKDVfXVqvpAklTVc6vqs0nSWns8yTuy\n/g77byb5xKTgDcCO07cBZrSjl52cSmvtVybc/3CS15xw/IUkL+xVF5O9/TNvzy2vvCXn7z1/6FKA\nAejbALNblE87YQTufuhuO10CAMxA+GbL9u/bn8M/FL4BALZL+GbLrth3hTPfAAAzEL7ZsisuuMKZ\nbwCAGQjfbNkV+67IkWNHhi4DAGC0hG+27KUXvzQ3XXfT0GUAAIyW8M2WLe9dzrWXXTt0GQAAoyV8\nAwBAJ8I3AAB0InwDAEAnwjcAAHQifDOVL333S3n3He8eugwAgFESvpnKM5aekYNHDg5dBgDAKAnf\nTGX/vv058siRtNaGLgUAYHSEb6ayvHc55+w5Jz947AdDlwIAMDrCN1Pbv29/7n/k/qHLAAAYHeGb\nqf380hMAAKazNHQBjM/Nr7g55+89f+gyAABGR/hmaiv7VoYuAQBglFx2AgAAnQjfAADQifANAACd\nCN8AANCJ8M22vO22t+XOI3cOXQYAwKgI32zLntqT+47eN3QZAACjInyzLXa5BACYnvDNtqzsW8mR\nY3a5BACYhvDNtthiHgBgesI327KyvJIjjxxJa23oUgAARkP4ZlsueOYFufcd9w5dBgDAqCwNXQDj\nVFW56LyLhi4DAGBUnPkGAIBOhG8AAOhE+AYAgE6Eb2bi004AALZO+GbbDv/wcK78wJVDlwEAMBrC\nN9t28S9dnMOPHM7jTzw+dCkAAKMgfLNte5f25sJzL8xD//vQ0KUAAIyC8M1MVpZXcv8j9w9dBgDA\nKAjfzGRl3/o28wAAnJ7wzUz2L+932QkAwBbZXp6ZvPe692bPWXuGLgMAYBSc+WYmgjcAwNYJ3wAA\n0InwDQAAnQjfAADQifDNzH76s5/mseOPDV0GAMDCE76Z2Ttvf2c+cs9Hhi4DAGDhCd/MzC6XAABb\nI3wzs/379uf+Y8I3AMDpCN/MbGXfivANALAFwjczc9kJAMDWCN/M7MJzL8y+Z+7ziScAAKexNHQB\njF9V5fC7Dg9dBgDAwnPmGwAAOhG+AQCgE+EbAAA6Eb4BAKAT4Zu5+NkTP8t9R+8bugwAgIUmfDMX\nj/7fo3nJrS9Ja23oUgAAFpbwzVws713O0llLOfqTo0OXAgCwsIRv5mZleSUPHHtg6DIAABaW8M3c\nXL58uW3mAQBOQfhmbi5fvtyZbwCAUxC+mZurL7o6Z+85e+gyAAAW1tLQBbB7vPHFbxy6BACAhebM\nNwAAdCJ8AwBAJ8I3AAB0InwDAEAnwjdzdejooXznR98ZugwAgIUkfDNXt37l1nz8Gx8fugwAgIUk\nfDNXNtoBAJhM+GauLl++PA/86IGhywAAWEjCN3PlzDcAwGTCN3N12fJlefDYg2mtDV0KAMDCEb6Z\nq/POOS83vPCG/Pj4j4cuBQBg4SwNXQC7zydf98mhSwAAWEjOfAMAQCfCNwAAdFK76Y1xVdV20+8D\nnDmqKq21GrqOnvRsYKxm6dnOfAMAQCfCNwAAdCJ8AwBAJ4N91GBV3ZLkhiQ/TXI4yVtaa/+zybgH\nkvwoyRNJjrfWrulZJwDr9G2A2Q155vuOJFe21q5KcijJX0wY90SS1dba1Ro4wKD0bYAZDRa+W2t3\nttae2Di8K8klE4ZWXB4DMDh9G2B2i9Ic35rk9gmPtSQHq+rLVfWHHWsCYDJ9G2AbdvSa76o6mOQ5\nJ96V9aZ8Y2vtMxtjbsz6NYEfm/A0L2utPVxVv5z1Zv6t1toXd7JugDOVvg2ws3Y0fLfWXnmqx6vq\nzUleneTlp3iOhze+/6CqPpXkmiQTm/iBAweevL26uprV1dVpSgboYm1tLWtra0OX8TS9+7aeDYzB\nPHv2YDtcVtX1Sd6X5NrW2tEJY85NclZr7dGqelbW3+zzV621OyaMt1saMEpj2OFy3n1bzwbGapae\nPWT4PpTknCQ/b+B3tdb+pKqem+SDrbXXVNVKkk9l/U+eS0k+2lq7+RTPqZEDozSS8D3Xvq1nA2M1\nyvC9EzRyYKzGEL7nTc8GxmqWnr0on3YCAAC7nvANAACdCN8AANCJ8A0AAJ0I3wAA0InwDQAAnQjf\nAADQifANAACdCN8AANCJ8A0AAJ0I3wAA0InwDQAAnQjfAADQifANAACdCN8AANCJ8A0AAJ0I3wAA\n0InwDQAAnQjfAADQifANAACdCN8AANCJ8A0AAJ0I3wAA0InwDQAAnQjfAADQifANAACdCN8AANCJ\n8A0AAJ0I3wAA0InwDQAAnQjfAADQifANAACdCN8AANCJ8A0AAJ0I3wAA0InwDQAAnQjfAADQifAN\nAACdCN8AANCJ8A0AAJ0I3wAA0InwDQAAnQjfAADQifANAACdCN8AANCJ8A0AAJ0I3wAA0InwDQAA\nnQjfAADQifANAACdCN8AANCJ8A0AAJ0I3wAA0InwDQAAnQjfAADQifANAACdCN8AANCJ8A0AAJ0I\n3wAA0InwDQAAnQjfAADQifANAACdCN8AANCJ8A0AAJ0I3wAA0InwDQAAnQjfAADQifANAACdCN8A\nANCJ8A0AAJ0I3wAA0InwDQAAnQjfAADQifANAACdCN8AANCJ8A0AAJ0I3wAA0InwDQAAnQjfAADQ\nifANAACdCN8AANCJ8A0AAJ0I3wAA0InwDQAAnQwWvqvqr6vqP6rqa1X1haq6aMK466vqv6rqvqr6\ns951ArBO3waY3ZBnvm9prb24tXZ1ks8luenkAVV1VpJ/TPK7Sa5M8oaqelHfMnevtbW1oUsYFfM1\nHfO1K+nbA/Kamo75mo756mew8N1ae/SEw2cleWKTYdckOdRae7C1djzJJ5K8tkd9ZwIvtOmYr+mY\nr91H3x6W19R0zNd0zFc/S0P+8Kr6myRvSnIsye9sMuTiJN894fh7WW/sAAxA3waYzY6e+a6qg1X1\n9RO+vrHx/YYkaa29p7V2aZKPJnnnTtYCwOnp2wA7q1prQ9eQqnpeks+31n79pPt/M8mB1tr1G8d/\nnqS11v52wvMM/8sAbFNrrYauYavm0bf1bGDMttuzB7vspKqe31r79sbh7yX51ibDvpzk+VV1WZKH\nk7w+yRsmPeeY/uECGJt59209GzgTDXnN981V9YKsv2HnwSR/lCRV9dwkH2ytvaa19nhVvSPJHVm/\nROZDrbXNmj0AO0/fBpjRQlx2AgAAZ4LR7XC5lc0bqurvq+pQVd1TVVf1rnGRnG6+quq6qjpWVV/d\n+HrPEHUuiqr6UFV9v6q+foox1teG082X9fVUVXVJVf1bVX1z442M75owbtesMT17evr21unZ09Gz\np7NjPbu1NpqvrP9n4dtJLktydpJ7krzopDGvSvK5jdu/keSuoete8Pm6LsltQ9e6KF9JfjvJVUm+\nPuFx62u6+bK+njofFyW5auP2eUnu3c09TM/esTnzuvrFXOjZ850va+up87EjPXtsZ763snnDa5P8\nU5K01u5Ocn5VPadvmQtjq5tdeNPThtbaF5M8cooh1tcJtjBfifX1pNbaf7fW7tm4/WjW37B48UnD\ndtMa07Onp29PQc+ejp49nZ3q2WML35tt3nDyJJw85qFNxpwptjJfSfJbG38q+VxV/Wqf0kbL+pqe\n9bWJqro862eg7j7pod20xvTs6enb82V9Tc/a2sQ8e/agO1yyEL6S5NLW2mNV9aok/5rkBQPXxO5h\nfW2iqs5L8i9J/rQ9dct22AqvK3aKtbWJeffssZ35fijJpSccX7Jx38ljnneaMWeK085Xa+3R1tpj\nG7dvT3J2VV3Qr8TRsb6mYH09XVUtZb2J/3Nr7dObDNlNa0zPnp6+PV/W1xSsrafbiZ49tvD95OYN\nVXVO1jdvuO2kMbcleVPy5E5rx1pr3+9b5sI47XydeF1SVV2T9Y+f/GHfMhdOZfI1b9bX002cL+tr\nUx9O8p+ttb+b8PhuWmN69vT07enp2dPRs6cz9549qstO2oTNG6rq7esPt1tba5+vqldX1beT/DjJ\nW4aseUhbma8kr6uqP05yPMlPkvz+cBUPr6o+lmQ1ybOr6jtJbkpyTqyvTZ1uvmJ9PUVVvSzJHyT5\nRlV9LUlL8pdZ/2SLXbfG9Ozp6dvT0bOno2dPZ6d6tk12AACgk7FddgIAAKMlfAMAQCfCNwAAdCJ8\nAwBAJ8I3AAB0InwDAEAnwjcAAHQifAMAQCf/DzQ49IL888OwAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x878d390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "tcr_min = 0.1\n",
    "tcr_max = 1.2\n",
    "tcr_slp = 2.\n",
    "tcr_off = -4.69e-1\n",
    "eps = 1e-8\n",
    "rhos = 2650.\n",
    "cb = 0.5\n",
    "bmz = np.arange(0.,2.5,0.05)\n",
    "zb = bmz/(rhos*cb)\n",
    "Nbed = len(bmz)\n",
    "tcr = np.zeros_like(bmz)\n",
    "tcr2 = np.zeros_like(bmz)\n",
    "tcr[0] = tcr_min\n",
    "tcr2[0] = tcr_min\n",
    "for k in range(1,Nbed) :\n",
    "    tcr[k] = tcr_min\n",
    "    if(bmz[k-1] >= eps) :\n",
    "        tcr[k] = np.exp((np.log(bmz[k-1])-tcr_off)/tcr_slp)\n",
    "        tcr2[k] = 0.4*bmz[k-1]**0.56\n",
    "    tcr[k] = min( max( tcr[k], tcr_min), tcr_max )\n",
    "    tcr2[k] = min( max( tcr2[k], tcr_min), tcr_max )\n",
    "    \n",
    "fig = plt.figure(figsize=(12,10))\n",
    "ax=fig.add_subplot(1,2,1)\n",
    "plt.plot(tcr,-bmz,tcr2,-bmz,'--')\n",
    "plt.axis([0, 2, -3, 0])\n",
    "plt.ylabel('Mass Depth (kg/m^2)')\n",
    "ax=fig.add_subplot(1,2,2)\n",
    "plt.plot(tcr,-zb*1000,tcr2,-zb*1000,'--')\n",
    "plt.axis([0, 2, -3, 0])\n",
    "plt.ylabel('Depth (mm)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Variable   Type       Data/Info\n",
      "-------------------------------\n",
      "Nbed       int        10\n",
      "bmz        ndarray    1: 1 elems, type `float64`, 8 bytes\n",
      "k          int        1\n",
      "tcr        ndarray    1: 1 elems, type `float64`, 8 bytes\n",
      "tcr_max    float      1.2\n",
      "tcr_min    float      0.1\n",
      "tcr_off    float      -0.469\n",
      "tcr_slp    float      2.0\n"
     ]
    }
   ],
   "source": [
    "whos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[  1800.   3000.   4200.   5400.   6600.   7800.   9000.  10200.] [ 0.01  0.08  0.16  0.24  0.32  0.4   0.48  0.6 ]\n",
      "10200.0\n",
      "(170L,)\n"
     ]
    }
   ],
   "source": [
    "# UGEMS stress time series\n",
    "dtt = np.array([ 30., 20., 20. ,20.,20.,20.,20.,20.])\n",
    "tt =  np.cumsum(60.*dtt)\n",
    "taulevels = np.array([0.01, 0.08, 0.16, 0.24, 0.32, 0.40, 0.48, 0.60])\n",
    "print tt,taulevels\n",
    "print(tt[-1])\n",
    "t = np.arange(0,tt[-1],60)\n",
    "#dt = np.diff(t)\n",
    "#dt = np.concatenate([dt,dt[:]])\n",
    "tau = np.zeros_like(t)\n",
    "for i in np.arange(len(tau)):\n",
    "    for j in np.arange(len(taulevels)):\n",
    "        if t[i]>tt[j]:\n",
    "            tau[i]=taulevels[j]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "E = np.zeros_like(t)\n",
    "m = np.zeros_like(t)\n",
    "M = 5e-4;\n",
    "for i in np.arange(len(t)-1):\n",
    "    # functional relationship of tauc(m) from lookup table\n",
    "    tc = np.interp(m[i], bmz, tcr)\n",
    "    E = 0.\n",
    "    if(tau[i]>=tc):\n",
    "        E = M*(tau[i]-tc)\n",
    "    dt = t[i+1]-t[i]\n",
    "    m[i+1]= m[i]+E*dt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x8dd45c0>]"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAs8AAAJPCAYAAACDwjlvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmUZGV9P/73MwuM7CSyyCqoQUFQAVFkGyUGXEGNIIiK\nJoB+xXiMEDQxih6jRkM0HvEXUUAFBZQdZR8YQQHZFXCAUZAdxbAzwzDMPL8/qkd7lp6url5u99zX\n65w6U1V9u+rTPPT0ez79ee4ttdYAAACDm9R0AQAAMFEIzwAA0CXhGQAAuiQ8AwBAl4RnAADokvAM\nAABd6io8l1L2LKXcWkq5vZRyxDI+vlsp5dFSyvV9t0+NfKkAANCsKYMdUEqZlOQbSXZPcn+Sa0op\nZ9Vab13i0MtqrW8dhRoBAGBc6KbzvEOS2bXWu2qt85OcnGSvZRxXRrQyAAAYZ7oJzxsmuaff43v7\nnlvSjqWUG0spPy2lbDki1QEAwDgy6NhGl65LskmtdU4p5Q1JzkzyNyP02gAAMC50E57vS7JJv8cb\n9T33Z7XWJ/vdP6+U8s1Syl/VWh/uf1wppQ6nWAAA6FatdcTHirsZ27gmyQtLKZuWUlZK8q4kZ/c/\noJSyXr/7OyQpSwbnRWqtbhP09pnPfKbxGtysXRtv1m9i36zfxL1Zu4l9Gy2Ddp5rrQtKKYcmuTCd\nsH1srXVWKeWQzofrMUn+vpTyoSTzk8xNsu+oVQwAAA3paua51np+ki2WeO5b/e4fneTokS0NAADG\nF1cYpGvTp09vugR6ZO0mNus3sVm/icvasSxlNGdClnqzUupYvh8AAO1USkltaMMgAAAQ4RkAALom\nPAMAQJeEZwAA6JLwDAAAXRKeAQCgS8IzAAB0SXgGAIAuCc8AANAl4RkAALokPAMAQJeEZwAA6JLw\nDAAAXRKeAQCgS8IzAAB0aUrTBQAAE99VVyWzZjVdxdh62cuSbbdtugrGmvAMAAzbYYcla66ZrLtu\n05WMnTXXFJ7bSHgGAIZtzpzk618XJlnxmXkGAIZtzpzkOc9pugoYfcIzADBsc+cmq6zSdBUw+oRn\nAGDY5swRnmkH4RkAGDZjG7SF8AwADEutnbEN4Zk2EJ4BgGGZNy+ZOjWZPLnpSmD0Cc8AwLDYLEib\nCM8AwLDYLEibCM8AwLDYLEibCM8AwLAY26BNhGcAYFh0nmkT4RkAGBadZ9pEeAYAhsWGQdpEeAYA\nhsXYBm0iPAMAw2JsgzYRngGAYdF5pk2EZwBgWHSeaRPhGQAYFhsGaRPhGQAYFmMbtInwDAAMi7EN\n2kR4BgCGReeZNhGeAYBh0XmmTYRnAGBYbBikTYRnAGBYjG3QJsIzADAsxjZoE+EZABgWYxu0ifAM\nAAzL3LnGNmgP4RkAGBadZ9pEeAYAhsWGQdpEeAYAhsWGQdpEeAYAhsXYBm0iPAMAPVuwIHnmmWTl\nlZuuBMaG8AwA9OzppzvzzqU0XQmMDeEZAOiZzYK0jfAMAPTMZkHaRngGAHpmsyBtIzwDAD0ztkHb\nCM8AQM+MbdA2wjMA0DOdZ9pGeAYAeqbzTNsIzwBAz2wYpG2EZwCgZ8Y2aBvhGQDombEN2kZ4BgB6\npvNM2wjPAEDPdJ5pG+EZAOiZDYO0jfAMAPTM2AZtIzwDAD0ztkHbCM8AQM90nmkb4RkA6JnOM20j\nPAMAPbNhkLYRngGAnhnboG2EZwCgZ8Y2aBvhGQDomc4zbSM8AwA903mmbYRnAKBnNgzSNsIzANAz\nYxu0jfAMAPTM2AZtIzwDAD2ZP7/z59SpzdYBY0l4BgB6YmSDNhKeAYCeGNmgjYRnAKAnOs+0kfAM\nAPRE55k2Ep4BgJ44xzNtJDwDAD0xtkEbCc8AQE+MbdBGwjMA0BOdZ9pIeAYAeqLzTBsJzwBAT2wY\npI2EZwCgJ8Y2aCPhGQDoibEN2kh4BgB6ovNMGwnPAEBPdJ5pI+EZAOiJDYO0kfAMAPTE2AZtJDwD\nAD0xtkEbCc8AQE90nmkj4RkA6InOM20kPAMAPbFhkDYSngGAnhjboI2EZwCgJ8Y2aCPhGQDoic4z\nbSQ8AwA90XmmjboKz6WUPUspt5ZSbi+lHLGc415ZSplfSnn7yJUIAIxHNgzSRoOG51LKpCTfSLJH\nkq2S7FdKefEAx30pyQUjXSQAML7U2uk8G9ugbbrpPO+QZHat9a5a6/wkJyfZaxnHfSTJqUn+OIL1\nAQDj0Lx5yUorJZMMgNIyU7o4ZsMk9/R7fG86gfrPSikbJNm71vraUspiHwOAFdXcuckjjzRdRTMe\nfVTXmXbqJjx342tJ+s9Cl4EOPPLII/98f/r06Zk+ffoIlQAAY2u//ZLLL09WXrnpSprxilc0XQH8\nxcyZMzNz5sxRf59Sa13+AaW8OsmRtdY9+x5/Ikmttf5nv2PuWHQ3yXOTPJXk4Frr2Uu8Vh3s/QBg\nothtt+Rzn+v8CYwvpZTUWgds6Paqm87zNUleWErZNMkDSd6VZL/+B9RaN190v5RyfJJzlgzOALCi\ncbYJaJ9Bw3OtdUEp5dAkF6azwfDYWuusUsohnQ/XY5b8lFGoEwDGHRcJgfYZdGxjRN/M2AYAK5DN\nN08uvrjzJzC+jNbYhhPMAECPdJ6hfYRnAOiRy1ND+wjPANAjGwahfYRnAOjB/PmdS1RPndp0JcBY\nEp4BoAdGNqCdhGcA6IHNgtBOwjMA9EDnGdpJeAaAHtgsCO0kPANAD4xtQDsJzwDQA2Mb0E7CMwD0\nQOcZ2kl4BoAemHmGdhKeAaAHxjagnYRnAOiBsQ1oJ+EZAHqg8wztJDwDQA90nqGdhGcA6IENg9BO\nwjMA9MDYBrST8AwAPTC2Ae0kPANAD3SeoZ2EZwDogc4ztJPwDAA9sGEQ2kl4BoAeGNuAdhKeAaAH\nxjagnYRnAOiBzjO0k/AMAD3QeYZ2Ep4BoAc2DEI7Cc8A0IO5c3WeoY2EZwDogc4ztJPwDABDVKvO\nM7SV8AwAQzRvXjJ1ajJ5ctOVAGNNeAaAITKyAe0lPAPAEBnZgPYSngFgiHSeob2EZwAYIlcXhPYS\nngFgiFxdENpLeAaAITK2Ae0lPAPAENkwCO0lPAPAEOk8Q3sJzwAwRDYMQnsJzwAwRDYMQnsJzwAw\nRMY2oL2EZwAYIhsGob2EZwAYIp1naC/hGQCGSHiG9hKeAWCIjG1AewnPADBEOs/QXsIzAAyRzjO0\nl/AMAEOk8wztJTwDwBAJz9BewjMADJGxDWgv4RkAhkjnGdpLeAaAIdJ5hvYSngFgiHSeob2EZwAY\nIuEZ2kt4BoAhMrYB7SU8A8AQLFiQPPNMMm1a05UATRCeAWAInn66E5xLaboSoAnCMwAMgXlnaDfh\nGQCGQHiGdhOeAWAIbBaEdhOeAWAIdJ6h3YRnABgCnWdoN+EZAIZA5xnaTXgGgCGYM0fnGdpMeAaA\nIZg7V+cZ2kx4BoAhMLYB7SY8A8AQ2DAI7SY8A8AQ6DxDuwnPADAENgxCuwnPADAENgxCuwnPADAE\nxjag3YRnABgCGwah3YRnABgCnWdoN+EZAIbAhkFoN+EZAIbAhkFoN+EZAIbA2Aa0m/AMAENgbAPa\nTXgGgCEwtgHtJjwDwBDoPEO7Cc8AMAQ6z9BuwjMADIENg9BuwjMADIGxDWg34RkAujR/fufPqVOb\nrQNozpSmCwBgYvn855Mf/ajpKpqxcGGyxhpJKU1XAjRFeAZgSK68Mjn44GTXXZuupBlrrtl0BUCT\nhGcAhmTu3GTLLZNttmm6EoCxZ+YZgCFxtgmgzYRnAIbE2SaANhOeARgSFwkB2kx4BmBIdJ6BNhOe\nARgSnWegzYRnAIbEhkGgzYRnALq2YEHyzDPJyis3XQlAM4RnALr29NOdeWdX2APaSngGoGs2CwJt\nJzwD0DWbBYG2E54B6JrNgkDbCc8AdG3uXGMbQLsJzwB0TecZaDvhGYCu2TAItJ3wDEDXbBgE2q6r\n8FxK2bOUcmsp5fZSyhHL+PhbSym/KqXcUEq5upSy08iXCkDTdJ6Btpsy2AGllElJvpFk9yT3J7mm\nlHJWrfXWfoddXGs9u+/4rZP8KMlLRqFeABqk8wy0XTed5x2SzK613lVrnZ/k5CR79T+g1jqn38PV\nkiwcuRIBGC9sGATarpvwvGGSe/o9vrfvucWUUvYupcxKck6SD4xMeQCMJ8Y2gLYbsQ2DtdYza60v\nSbJ3ks+P1OsCMH4Y2wDabtCZ5yT3Jdmk3+ON+p5bplrrz0spm5dS/qrW+vCSHz/yyCP/fH/69OmZ\nPn1618UC0Kw5c5I112y6CoClzZw5MzNnzhz19ym11uUfUMrkJLels2HwgSRXJ9mv1jqr3zEvqLX+\nru/+tknOqrVuvIzXqoO9HwDj10c/mmy+eedPgPGslJJaaxnp1x2081xrXVBKOTTJhemMeRxba51V\nSjmk8+F6TJJ3lFLem+SZJHOT7DPShQLQPBsGgbbrZmwjtdbzk2yxxHPf6nf/y0m+PLKlATDe2DAI\ntJ0rDALQNRsGgbYTngHoms4z0HbCMwBd03kG2k54BqBrNgwCbSc8A9A1YxtA2wnPAHTN2AbQdsIz\nAF3TeQbaTngGoGs6z0DbCc8AdM2GQaDthGcAujJ/flJrMnVq05UANEd4BqArRjYAhGcAumSzIIDw\nDECXdJ4BhGcAumSzIIDwDECXjG0ACM8AdMnYBoDwDECXdJ4BhGcAuqTzDCA8A9AlGwYBhGcAumRs\nA0B4BqBLxjYAhGcAuqTzDCA8A9AlnWcA4RmALtkwCCA8A9CluXONbQAIzwB0RecZQHgGoEs2DAII\nzwB0yYZBAOEZgC7pPAMIzwB0SecZQHgGoEs2DAIIzwB0ydgGgPAMQJeMbQAIzwB0SecZQHgGoEs6\nzwDCMwBdqFXnGSARngHowrx5ydSpyeTJTVcC0CzhGYBBGdkA6BCeARiUkQ2ADuEZgEHpPAN0CM8A\nDMrVBQE6hGcABmVsA6BDeAZgUMY2ADqEZwAGpfMM0CE8AzAonWeADuEZgEHZMAjQITwDMChjGwAd\nwjMAgzK2AdAhPAMwKJ1ngA7hGYBB6TwDdAjPAAzKhkGADuEZgEEZ2wDoEJ4BGJSxDYAO4RmAQek8\nA3QIzwAMSucZoEN4BmBQNgwCdAjPAAzK2AZAh/AMwKCMbQB0CM8ADErnGaBDeAZgUDrPAB3CMwCD\n0nkG6BCeARiUzjNAx5SmCwCYCC66qHNrq3nzkmnTmq4CoHnCM0AXjj02mTw5ednLmq6kGd/5TlJK\n01UANE94BujCnDnJP/xDstdeTVcCQJPMPAN0wcwvAInwDNAVZ5sAIBGeAbqi8wxAIjwDdGXOHOEZ\nAOEZoCvGNgBIhGeArhjbACARngG6ovMMQCI8Awyq1k7nWXgGQHgGGMS8eclKK3WuMAhAuwnPAIMw\nsgHAIsIzwCBsFgRgEeEZYBA6zwAsIjwDDELnGYBFhGeAQbi6IACLCM8AgzC2AcAiwjPAIIxtALCI\n8AwwCJ1nABYRngEGofMMwCLCM8AgbBgEYBHhGWAQxjYAWER4BhiEsQ0AFhGeAQah8wzAIsIzwCB0\nngFYRHgGGIQNgwAsIjwDDMLYBgCLCM8AgzC2AcAiwjPAIHSeAVhEeAYYhM4zAIsIzwCDsGEQgEWE\nZ4BBGNsAYBHhGWAQxjYAWER4BhiEzjMAiwjPAIPQeQZgEeEZYBA6zwAsIjwDLMfChckzzyTTpjVd\nCQDjgfAMsBxz53a6zqU0XQkA44HwDLAcRjYA6E94BlgOmwUB6K+r8FxK2bOUcmsp5fZSyhHL+Pj+\npZRf9d1+XkrZeuRLBRh7Os8A9DdoeC6lTEryjSR7JNkqyX6llBcvcdgdSXattb4syeeTfHukCwVo\ngs4zAP1103neIcnsWutdtdb5SU5Oslf/A2qtV9VaH+t7eFWSDUe2TIBmzJkjPAPwF92E5w2T3NPv\n8b1Zfjj+xyTnDacogPHC2AYA/U0ZyRcrpbw2yfuT7DySrwvQFGMbAPTXTXi+L8km/R5v1PfcYkop\n2yQ5JsmetdZHBnqxI4888s/3p0+fnunTp3dZKsDY03kGmBhmzpyZmTNnjvr7lFrr8g8oZXKS25Ls\nnuSBJFcn2a/WOqvfMZskmZHkPbXWq5bzWnWw9wMYT44/Prnsss6fAEwcpZTUWkf8EleDdp5rrQtK\nKYcmuTCdGelja62zSimHdD5cj0ny70n+Ksk3Syklyfxa6w4jXSzAWLNhEID+upp5rrWen2SLJZ77\nVr/7ByU5aGRLA2iesQ0A+nOFQYDlsGEQgP6EZ4Dl0HkGoD/hGWA5dJ4B6E94BlgOGwYB6E94BlgO\nYxsA9Cc8AyyHsQ0A+hOeAZZD5xmA/oRngOXQeQagP+EZYDlsGASgP+EZYDmMbQDQn/AMsBzGNgDo\nT3gGWA6dZwD6E54BlkPnGYD+hGeA5bBhEID+hGeAAcyfn9SaTJ3adCUAjBfCM8AAjGwAsCThGWAA\nNgsCsCThGWAAOs8ALEl4BhiAzYIALEl4BhiAsQ0AliQ8AwzA2AYASxKeAQag8wzAkoRngAHoPAOw\nJOEZYAA2DAJLqrXmZ7//Wd57xnvzP1f9T9Pl0IApTRcAMF7NnWtsA+i4/4n7870bv5fjbjwuK01e\nKf/win/I/lvv33RZNEB4BhiAzjO027MLn81Pb/9pjr3h2Fx+9+X5+5f8fU542wl51YavSiml6fJo\niPAMMAAbBqGdHp77cL5z/Xdy9DVHZ8PVN8xB2x6UH77jh1ltpdWaLo1xQHgGGIANg9Aut/zxlnz9\nl1/Pj37zo7zlb96S0/c5PdttsF3TZTHOCM8AA5gzJ1lnnaarAEZTrTXnzj43X/vl13LzH2/OB7f7\nYGZ9eFbWX239pktjnBKeAQag8wwrrgULF+S0WaflC5d/IUny8R0/nn222icrT1m54coY74RngAHY\nMAgrnvkL5ueHN/0wX/z5F7PWtLXy+dd9Pm960ZtsAKRrwjPAAGwYhBXHvGfn5bs3fjdf+sWXstla\nm+XoNx6d1232OqGZIROeAQZgbAMmvgULF+TEX5+YT8/8dLZcZ8uc+LYTs9MmOzVdFhOY8AwwAJ1n\nmLhqrTnvt+flExd/IquvvHp+8PYfZOdNdm66LFYAwjPAAHSeYWK6+r6rc8TFR+TBJx/Ml3b/Ut66\nxVuNZzBihGeAAdgwCBPLHY/ckSMuPiJX3nNljpx+ZA58+YGZMknUYWRNaroAgPHK2AZMDE8/+3Q+\n97PPZYdv75CXr/fy3P6R2/OP2/6j4Myo8H8VwACMbcD4d+7sc/NP5/1TXr7+y3P9IddnkzU3abok\nVnDCM8AAdJ5h/Pr9o7/Pxy74WG7+4835xhu/kT1fuGfTJdESwjMwoAcfTG65pekqmvPkkzrPMN7M\nXzA/X7niKznqyqPysVd/LCe946RMmzKt6bJoEeEZGNCRRyY/+1mywQZNV9KMN7xBeIbx5KY/3JT3\nnfm+rLfaern2oGuz2dqbNV0SLSQ8AwN68snk3/4tOeCApisB2uzZhc/mK7/4Sv77qv/Ol3b/Uj7w\nig849RyNEZ6BAZn5BZo266FZOfCsA7PGymvkuoOvsyGQxjlVHTAgZ5sAmrJg4YL81xX/lV2/u2ve\n//L358IDLhScGRd0noEB6TwDTbjv8fuy/+n7Z1KZlKv/8WqzzYwrOs/AgHSegbF2wW8vyPbf3j57\nvGCPzHjvDMGZcUfnGRiQy1MDY+XZhc/myJlH5rs3fjcnv+Pk7Pb83ZouCZZJeAYGZGwDGAsPPPFA\n9jttv6w0eaVcf8j1WXfVdZsuCQZkbAMYkLENYLRdfMfF2e6Y7bL7ZrvnvHefJzgz7uk8AwPSeQZG\nS601X7j8C/nmtd/MD97+g7x2s9c2XRJ0RXgGBqTzDIyGOfPn5ANnfSB3PnpnrjnommyweksvY8qE\nZGwDWKb585OFC5OpU5uuBFiR3PPYPdnl+F0ydfLU/OzAnwnOTDjCM7BMc+d2RjZcARcYKTc+eGN2\nPHbH7LvVvvn+3t/PtCnTmi4JhszYBrBMRjaAkXTR7y7Ku09/d45+49F551bvbLoc6JnwDCyTzYLA\nSDnhVyfksIsOy2n7nJZdNt2l6XJgWIRnYJl0noGR8LWrvpavXvXVXPq+S7PlOls2XQ4Mm/AMLJOr\nCwLDsehUdN/91Xdz2YGXZdO1Nm26JBgRwjOwTMY2gF7VWvPJGZ/MT27/SS478LI8b/XnNV0SjBjh\nGVgmYxtAL2qt+ecL/jmX3X1ZZh44M89d5blNlwQjSngGlsnYBjBUtdYccfER+dldP8uM987I2s9Z\nu+mSYMQJz8AyLTrPM0C3Pn3pp3P+b8/Ppe+7VHBmhSU8A8uk8wx04+lnn841912Tk28+OTPvmplL\n33dp/nqVv266LBg1wjOwTDYMAsvy7MJnc/V9V+fiOy7OjDtn5Lr7r8tL1nlJdt1k18x474ysu+q6\nTZcIo0p4BpbJhkEg6cwxz354di763UW56I6LMvP3M/P8tZ6f12/++nxy509mp413yuorr950mTBm\nhGdgmYxtQHvNXzA/l999ec669aycffvZeXbhs3n95q/Pvlvtm2+9+VtZb7X1mi4RGiM8A8s0Z06y\nzjpNVwGMlSefeTIX/PaCnHnbmTl39rl5wdovyF5b7JVz9jsnW62zVUopTZcI44LwDCyTsQ1Y8T39\n7NP56e0/zYk3nZgZd8zIjhvvmL222Ctf3P2L2WiNjZouD8Yl4RlYJhsGYcW0sC7M5XddnhN/fWJO\nm3VaXvG8V+SArQ/I8Xsdn7WmrdV0eTDuCc/AMuk8w4rldw//LsfecGx+cNMPsubKa+Y927wnv/7Q\nr3WYYYiEZ2CZbBiEiW/BwgU5d/a5+ea138y191+b973sfTlnv3OyzXrbNF0aTFjCM7BMxjZg4nro\nqYdy7A3H5n+v/d+st9p6+X/b/7+cvs/pec5U39QwXMIzsEzGNmDiueWPt+TLV3w5Z916Vt7+krfn\n1H1OzfYbbN90WbBCEZ6BZdJ5honj+geuz39c/h/5+d0/z0df9dH87p9+5xLZMEqEZ2CZdJ5h/PvF\n3b/I5y//fG76w0057DWH5ft7fz+rrrRq02XBCk14BpbJhkEYvy6989J89mefzd2P3Z1P7PyJnLnv\nmVl5yspNlwWtIDwDy2RsA8afGx64IZ+c8cnMfnh2PrPbZ7L/1vtnyiQ/ymEs+Y4DlsnYBowfdz92\ndz4545O55M5L8qldPpWDtjsoK01eqemyoJUmNV0AMP7UqvMM48Gc+XPy2ZmfzSu+9Yq8YO0XZPZH\nZufDO3xYcIYG6TwDS3nmmWTq1GTy5KYrgXaqteaMW8/Ixy74WF690atz/cHXZ9O1Nm26LCDCM7AM\nNgtCc37/6O9z6LmH5o5H7sj39/5+dnv+bk2XBPRjbANYipENGHsLFi7If13xX9n+mO3zmo1fkxs/\neKPgDOOQzjOwFJsFYWzd+qdb8/6z3p+VJ6+cq/7xqrzwr17YdEnAAHSegaXoPMPYWFgX5qgrjsrO\nx+2cA7Y+IJe87xLBGcY5nWdgKTrPMPrue/y+vPfM9+aZBc/k6oOuzuZrb950SUAXdJ6BpdgwCKPr\nrFvPyrbHbJvpm07PzPfNFJxhAtF5BpZibANGxzMLnsnhFx6es28/O2fse0Zes/Frmi4JGCLhGViK\nsQ0Yefc+fm/2+fE+WWfVdXLDITdkrWlrNV0S0ANjG8BSdJ5hZF1y5yV55bdfmb222Ctn7HuG4AwT\nmM4zsBSdZxg5R199dD5/+efzg7f/IK/b7HVNlwMMk/AMLMWGQRi++Qvm56PnfzSX3XVZrvjAFdls\n7c2aLgkYAcIzsBRjGzA8jz79aN7xo3fkOVOekyv+4YqssfIaTZcEjBAzz8BSjG1A7+59/N7scvwu\neek6L81Z7zpLcIYVjPAMLMXYBvTmlj/ekp2O2ynv3ea9+dqeX8vkSZObLgkYYcY2gKXMnWtsA4bq\ninuuyNtOeVuO+rujcsA2BzRdDjBKhGdgKTrPMDSX3HlJ9j1135z4thOzxwv3aLocYBQJz8BSbBiE\n7p07+9wceOaBOfWdp2a35+/WdDnAKBOegaXYMAjLN3/B/Mx+eHYuvfPSfO6yz+Xs/c7Oqzd6ddNl\nAWNAeAaWYmwD/uLxeY/nVw/+Ktc/cH1uePCG3PDgDbn9/27PxmtsnK3X2zrnv/v8vOJ5r2i6TGCM\nCM/AUoxt0FZz58/NdQ9clyvvuTLXPXBdrn/g+tz3xH3Zet2t84r1X5GdNt4ph+5waLZcZ8usMtW/\nMKGNhGdgKcY2aIt7HrsnV957Za6454pcee+VufmPN2fLdbbMjhvtmDe96E351K6fyouf++JMmeTH\nJdDhbwNgKTrPrKgeeOKBXHLnJZ3b7y/JU888ldds/JrsuNGO+crrv5LtN9heRxlYLuEZWIrOMyuK\nx+c9nhl3zMiMO2fkkjsvyYNPPpjpz5+e3TfbPYe95rC8+LkvTiml6TKBCUR4BpZiwyAT2Z2P3Jlz\nbj8n59x+Tq6696rsuNGOef3mr88JbzshL1//5a76BwxLqbUOflApeyb5WjqX8z621vqfS3x8iyTH\nJ9k2yb/WWv97gNep3bwf0Kxp05JHHjG6wcSwsC7ML+/9Zc667az85Paf5KE5D+VNL3pT3vw3b87r\nN399Vl959aZLBBpQSkmtdcR/tTRoeC6lTEpye5Ldk9yf5Jok76q13trvmOcm2TTJ3kkeEZ5h4lq4\nMJkyJVmwIPHbbMarWmt+9Ydf5aSbTsrJt5ycVaeumre9+G15yxZvyQ4b7pBJZVLTJQING63w3M3Y\nxg5JZtda7+or5OQkeyX5c3iutf4pyZ9KKW8e6QKBsTV3bqfzLDgzHs3+v9k56eaTctLNJ2Xu/LnZ\n76X75Zz9zsnW625tdhkYE92E5w2T3NPv8b3pBGpgBWSzIOPNk888mVNuPiXfvv7b+f2jv88+W+2T\n4956XF690asFZmDMjfmGwSOPPPLP96dPn57p06ePdQnActgsyHgx66FZ+fovv55Tbjklu2y6Sz61\n66ey5wvwbAjhAAASf0lEQVT3dM5lYJlmzpyZmTNnjvr7dDPz/OokR9Za9+x7/IkkdclNg30f+0yS\nJ8w8w8R1663JXnslt93WdCW0Ua01M+6cka9e9dVcd/91+eD2H8xB2x6UDdfYsOnSgAmmyZnna5K8\nsJSyaZIHkrwryX7LOd7v0GACM7ZBE2qtOe+35+XImUfmyWeezMd3/HhO2+e0TJsyrenSABYzaHiu\ntS4opRya5ML85VR1s0oph3Q+XI8ppayX5NokqydZWEr5aJIta61PjmbxwMhzdUHG2oW/uzCfuuRT\nmfvs3Hxmt8/k7S95u7NlAONWV4Njtdbzk2yxxHPf6nf/D0k2HtnSgCboPDNWbv7jzTnswsNy56N3\n5vOv/XzeseU7hGZg3LPrAliMDYOMtofnPpx/nfGvOX3W6fn3Xf89H9z+g5k6eWrTZQF0xT/xgcUY\n22C01Frzw5t+mK2+uVWmTJqS2w69LR951UcEZ2BC0XkGFmNsg9Fw92N35+BzDs79T9yfM/c9M6/a\n6FVNlwTQE51nYDE6z4y0k246Kdsfs3123XTXXHfwdYIzMKHpPAOL0XlmpDz29GP58LkfzrX3X5vz\nDzg/2z5v26ZLAhg2nWdgMTYMMhJ+/YdfZ7tjtsvqK62e6w+5XnAGVhg6z8Bi5sxJ1lyz6SqYyE65\n+ZQcet6h+doeX8u7t3l30+UAjCjhGVjM3LnJ857XdBVMRAsWLsgnZ3wyP/7Nj3PRey7Ky9d/edMl\nAYw44RlYjA2D9OLpZ5/Oe854Tx566qFce9C1+etV/rrpkgBGhZlnYDE2DDJUj8x9JH93wt9lUpmU\nCw64QHAGVmjCM7AYGwYZivsevy+7HL9LtnvedjnpHSdl5SkrN10SwKgSnoHFGNugW/c+fm+mf296\n3rPNe/LVPb+aScWPFGDFZ+YZWIyxDbpx7+P35rXfe20O3vbgHL7T4U2XAzBmtAmAxRjbYDCCM9Bm\nOs/AYubONbbB4v4050+Z+fuZ+fndP88V91yRWx66JZ+d/tkc9prDmi4NYMwJz8BidJ6Z9+y8XHHP\nFbnwdxfmojsuyuyHZ2eXTXbJLpvskqP+7qhsv8H2ec5U/8IC2kl4BhZjw2A73f3Y3Tnz1jNz/m/P\nz8/v/nm2XGfLvH7z1+ere3w1r9roVVlp8kpNlwgwLpRa69i9WSl1LN8PevGFLySf+UzTVTRn6tTk\noYeSVVdtuhJG212P3pWTbz45p806LXc8ckfeusVb86YXvSmv2+x1Wfs5azddHsCwlFJSay0j/rrC\nMyzukEOSbbZJDj646UqaMWlSMnly01UwWh57+rH8+Dc/zgm/PiG3/PGW/P2Wf593bvnO7Pb83TJl\nkl9GAiuO0QrP/qaEJTz1VLLGGp0OLKworrv/uvzvtf+bU2edmtdt9rp87NUfyxtf9EbjGABDJDzD\nEp56ysgCK4YFCxfk9Fmn56grj8qDTz6Yg7c7OLM+PCvrr7Z+06UBTFjCMyxBeGaimzt/bo674bgc\ndeVRed7qz8sndv5E3vI3b8nkSeZxAIZLeIYlCM9MVPOenZdjbzg2X7j8C9n2edvmxLefmNds/Jqm\nywJYoQjPsAThmYlmYV2YE351Qj4989PZap2tcua7zsz2G2zfdFkAKyThGZYgPDORXHnPlfno+R/N\npDIpP3z7D7PTJjs1XRLACk14hiUIz0wEDz31UP75wn/OpXdemi/97Zey/9b7Z1KZ1HRZACs8f9PC\nEoRnxrNaa37w6x9k6/9v66y36nq59dBbc8A2BwjOAGNE5xn6qVV4Zvx64IkHctA5B+Xux+7OOfud\nk1du+MqmSwJoHa0K6OeZZzpX2HOBFMabn97+02x7zLbZfoPtc+3B1wrOAA3ReYZ+dJ0Zb+Y9Oy9H\nXHxEzrj1jPz4nT/Ozpvs3HRJAK0mPEM/wjPjyX2P35e3nfK2bLzmxrnxkBuz9nPWbrokgNYztgH9\nCM+MF1fde1V2+M4OeduL35ZT33mq4AwwTug8Qz/CM+PB9278Xg6/6PAct9dxefPfvLnpcgDoR3iG\nfoRnmlRrzRd//sV8+/pvZ+aBM7PlOls2XRIASxCeoR/hmaYsrAvz8Qs+novvvDi/+MAvssHqGzRd\nEgDLIDxDP8IzTZi/YH4+cPYHcucjd+ayAy8z3wwwjgnP0M+cOcIzY+vZhc/mgDMOyBPznsiF77kw\nq0xdpemSAFgO4Rn60XlmLC1YuCDvP+v9efTpR3PWu87KtCnTmi4JgEEIz9CP8MxYWVgX5uBzDs59\nj9+Xn+z/E8EZYIIQnqEf4Zmx8s8X/HNu+7/bcv4B5xvVAJhAhGfo56mnkuc+t+kqWNH9z1X/k4vu\nuCi/+MAvstpKqzVdDgBD4AqD0I/OM6PtzFvPzJev+HLO3f/crDVtrabLAWCIdJ6hH+GZ0XTNfdfk\noHMOynnvPi+brrVp0+UA0AOdZ+hHeGa0PPDEA9n7lL3znbd8J9tvsH3T5QDQI+EZ+hGeGQ3zF8zP\nvqfum0O2OyR7vXivpssBYBiEZ+hHeGY0/OuMf81qK62WT+36qaZLAWCYzDxDP8IzI+30Wafnx7/5\nca47+LpMKvoVABOd8Az9CM+MpN8+/Nt88CcfzE/3/2n+epW/brocAEaANgj089RTySquV8EImL9g\nft59+rvz77v+e1654SubLgeAESI8Qz86z4yUz/7ss3nuKs/NoTsc2nQpAIwgYxvQj/DMSLjsrsty\n3A3H5cYP3phSStPlADCCdJ6hz4IFyfz5ybRpTVfCRPbI3EfynjPek2PfemzWXXXdpssBYIQJz9Bn\n0byzRiG9qrXmgz/9YPbaYq+84UVvaLocAEaBsQ3oY2SD4frer76X3zz0m3xv7+81XQoAo0R4hj7C\nM8Px24d/m8MvOjyXvu/STJti9gdgRWVsA/oIz/Rq0WnpPr3rp/PSdV/adDkAjCLhGfoIz/TKaekA\n2sPYBvQRnumF09IBtIvOM/QRnhkqp6UDaB/hGfoIzwyF09IBtJOxDegjPDMUTksH0E7CM/QRnumW\n09IBtJexDegjPNMNp6UDaDfhGfoIz3TDaekA2s3YBvR56qlk/fWbroLxzGnpANB5hj46zyyP09IB\nkAjP8GfCMwOZ9dCsvOWktzgtHQDCMywiPLOkJ+Y9kX+56F+y63d3zT5b7ZP/3uO/my4JgIaZeYY+\nwjOL1Frzo1t+lMMuOiy7b7Z7bv7QzVlvtfWaLguAcUB4hj7CM0nym4d+k4+c95H8ac6fctI7TsrO\nm+zcdEkAjCPGNqCP8Nxui0Y0dvvubtl7i71z3cHXCc4ALEXnGfoIz+1kRAOAoRCeoY/w3D5GNAAY\nKmMb0Ed4bg8jGgD0SucZktSazJ2brLJK05UwmoxoADBcwjOkE5ynTUsm+V3MCsuIBgAjQVSAGNlY\nkRnRAGAk6TxDhOcVkRENAEaD8AwRnlc0RjQAGC3GNiDC84rCiAYAo03nGSI8T3RGNAAYK8IzRHie\nyIxoADCWjG1AhOeJyIgGAE3QeYYIzxOJEQ0AmiQ8Q4TnicKIBgBNM7YBEZ7HOyMaAIwXOs+QTnhe\nZZWmq2BJRjQAGG+EZ0gnPG+4YdNV0N+iEY3/m/N/OfkdJ2enTXZquiQAMLYBibGN8eSJeU/k8AsP\n//OIxrUHXys4AzBuCM8Q4Xk8qLXmlJtPyUuOfkkemvNQbv7QzfnIqz6SKZP8ggyA8cNPJYjw3LT+\nIxqn/P0pOs0AjFs6zxDhuSlGNACYaIRniPA81oxoADBR+UkFEZ7HkhENACYynWeI8DwWjGgAsCLQ\neSZJcvzxyX/8R9NVNOeuu5K11mq6ihXTogudfPzCj+dvN/9bFzoBYEIrtdaxe7NS6li+H93bfvvk\n8MOT7bZrupJmTJuWbLRR01WsePqPaBz9xqN1mgEYM6WU1FrLiL+u8Mz99ycvfWnyhz8kU6c2XQ0r\ngifmPZHP/exz+e6vvptP7/rpfOiVH7IZEIAxNVrh2cwz+elPkz32EJwZPmfRAGBF5ycaOeec5F3v\naroKJjpn0QCgDYxttNzcucl663U2zK29dtPVMBEZ0QBgPDK2waiYMSPZdlvBmaEzogFAG3X1U66U\nsmeSr6UTto+ttf7nMo75epI3JHkqyYG11htHslBGx09+krzlLU1XwURjRAOAthq081xKmZTkG0n2\nSLJVkv1KKS9e4pg3JHlBrfVFSQ5J8r+jUCsjrNZOeH7zm7s7fubMmaNaD6NnpNbOhU6a4XtvYrN+\nE5e1Y1m6GdvYIcnsWutdtdb5SU5OstcSx+yV5PtJUmv9ZZI1SymugjDO3XBDssoqyRZbdHe8v0Qm\nruGunRGNZvnem9is38Rl7ViWbn7ybZjknn6P700nUC/vmPv6nvvDsKpjVBnZoBtGNADgL8a8bbTe\nx6S18eLRR5Pttk3eclJ3x99202257qTrRrcoRkWvazd/wfxc98B1zqIBAH0GPVVdKeXVSY6ste7Z\n9/gTSWr/TYOllP9Ncmmt9ZS+x7cm2a3W+oclXst56gAAGBOjcaq6btpI1yR5YSll0yQPJHlXkv2W\nOObsJB9Ockpf2H50yeCcjM4XAAAAY2XQ8FxrXVBKOTTJhfnLqepmlVIO6Xy4HlNrPbeU8sZSym/T\nOVXd+0e3bAAAGHtjeoVBAACYyMbsCoOllD1LKbeWUm4vpRwxVu/LwEopG5VSLiml3FJKuamU8k99\nz69dSrmwlHJbKeWCUsqa/T7nk6WU2aWUWaWUv+v3/LallF/3re/Xmvh62qiUMqmUcn0p5ey+x9Zu\ngiilrFlK+XHfetxSSnmV9ZsYSikfK6Xc3Pff/QellJWs3fhVSjm2lPKHUsqv+z03YuvVt/4n933O\nlaWUTcbuq1vxDbB+X+5bnxtLKaeVUtbo97HRX79a66jf0gnpv02yaZKpSW5M8uKxeG+35a7L+kle\n3nd/tSS3JXlxkv9M8i99zx+R5Et997dMckM64z7P71vTRb+9+GWSV/bdPzfJHk1/fW24JflYkhOT\nnN332NpNkFuS7yZ5f9/9KUnWtH7j/5ZkgyR3JFmp7/EpSd5n7cbvLcnOSV6e5Nf9nhux9UryoSTf\n7Lu/b5KTm/6aV6TbAOv3t0km9d3/UpIvjuX6jVXnuZsLrTDGaq0P1r7LqNdan0wyK8lG6azN9/oO\n+16SvfvuvzWd/6merbX+PsnsJDuUUtZPsnqt9Zq+477f73MYJaWUjZK8Mcl3+j1t7SaAvi7JLrXW\n45Okb10ei/WbKCYnWbWUMiXJc9K5toG1G6dqrT9P8sgST4/kevV/rVOT7D7iX0SLLWv9aq0X11oX\n9j28Kp3skozR+o1VeF7WhVY2HKP3pgullOen8y+7q5KsV/vOllJrfTDJun2HDXQxnA3TWdNFrO/Y\n+GqSw5P037hg7SaGzZL8qZRyfN/YzTGllFVi/ca9Wuv9SY5Kcnc66/BYrfXiWLuJZt0RXK8/f06t\ndUGSR0spfzV6pbOED6TTSU7GaP3GbOaZ8auUslo6/9r6aF8HesldpHaVjjOllDcl+UPfbw6WdwpI\nazc+TUmybZKja63bpnOWok/E9964V0pZK51O1abpjHCsWkp5d6zdRDeS6+W0vGOklPJvSebXWru8\n3Ft3LzvYAWMVnu9L0n8Ae6O+52hY368dT01yQq31rL6n/1BKWa/v4+sn+WPf8/cl2bjfpy9ax4Ge\nZ/TslOStpZQ7kpyU5HWllBOSPGjtJoR7k9xTa7227/Fp6YRp33vj398muaPW+nBfl+qMJK+JtZto\nRnK9/vyxUsrkJGvUWh8evdJJklLKgemMLu7f7+kxWb+xCs9/vtBKKWWldC60cvYYvTfLd1yS39Ra\n/6ffc2cnObDv/vuSnNXv+Xf17UzdLMkLk1zd9yuvx0opO5RSSpL39vscRkGt9V9rrZvUWjdP5/vp\nklrre5KcE2s37vX9uvieUsrf9D21e5Jb4ntvIrg7yatLKdP6/pvvnuQ3sXbjXcniHcWRXK+z+14j\nSd6Z5JJR+yraa7H1K6Xsmc7Y4ltrrfP6HTc26zeGuyX3TOdsDrOTfGKs3tdtuWuyU5IF6Zz95IYk\n1/et018lubhvvS5Msla/z/lkOrtXZyX5u37Pb5fkpr71/Z+mv7Y23ZLslr+cbcPaTZBbkpel01i4\nMcnp6Zxtw/pNgFuSz/Stw6/T2Wg01dqN31uSHya5P8m8dP7x8/4ka4/UeiVZOcmP+p6/Ksnzm/6a\nV6TbAOs3O8ldfbnl+vSdLWOs1s9FUgAAoEs2DAIAQJeEZwAA6JLwDAAAXRKeAQCgS8IzAAB0SXgG\nAIAuCc8AANAl4RkAALr0/wM+5G1QtA09cAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x8be5518>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(12,10))\n",
    "plt.plot(t,tau)\n",
    "plt.plot(t,m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x903ecf8>]"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAs4AAAJPCAYAAABsAFJRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH6xJREFUeJzt3X2spnV95/HPlxksUJDGrCGRB2utBSVtrdlSWtt4GroR\ncQtkt7vr0KzRzTZkK0pK29A0rE6T3T82rVtrcFfIogmNDW7dprpZa9ime0isqShKa+3w4NpQQMoW\nrVXAh2H47h/ngMfDYc53zpw598z4eiUn53r4Xdf9m+HK4c3FdZ+7ujsAAMDBnbDoCQAAwLFAOAMA\nwIBwBgCAAeEMAAADwhkAAAaEMwAADIzCuaourqq7quqeqrp2g/3nVtXHqurrVXXNun2/VFV/WVV/\nUVXvq6rnbNfkAQBgp2wazlV1QpLrk7w6yflJ9lTVeeuGfTHJm5P85rpjX7C6/RXd/UNJdid53TbM\nGwAAdtTkjvMFSe7t7vu6e3+SW5JctnZAdz/S3XckeWKD43cl+e6q2p3klCRfOMw5AwDAjpuE85lJ\n7l+z/sDqtk119xeSvD3J3yR5MMmXu/uPD3WSAACwaEf0zYFV9T1ZuTv9wiQvSHJqVV1xJF8TAACO\nhN2DMQ8mOWfN+lmr2yZ+Jsnnu/tLSVJVf5DkJ5L83vqBVdXDcwIAwJZ1d23luEk4fyLJ91fVC5M8\nlJU39+05yPi1E/mbJBdW1UlJvpHkotXzbahbO/Mte/fuzd69exc9DY4yrgs24rpgI64LNlK1pWZO\nMgjn7j5QVVcluTUrj3bc1N37qurKld19Y1WdkeSTSU5L8mRVXZ3kZd19e1V9IMmnk+xf/X7jlmcL\nAAALMrnjnO7+SJJz1227Yc3yw0nOfpZjfyPJbxzGHAEAYOF8ciBHraWlpUVPgaOQ64KNuC7YiOuC\n7VZHy3PFVdVHy1wAADg+VdWW3xzojjMAAAwIZwAAGBDOAAAwIJwBAGBAOAMAwIBwBgCAAeEMAAAD\nwhkAAAaEMwAADAhnAAAW5ktfSn7lVxY9ixnhDADAwjz+eHLLLYuexYxwBgCAAeEMAAADwhkAAAaE\nMwAADAhnAAAYEM4AADAgnAEAYEA4AwDAgHAGAIAB4QwAAAPCGQAABoQzAAAMCGcAABgQzgAAMCCc\nAQBgQDgDAMCAcAYAgAHhDAAAA8IZAAAGhDMAAAwIZwAAGBDOAAAwIJwBAGBAOAMAwIBwBgCAAeEM\nAAADwhkAAAaEMwAADAhnAAAYEM4AADAgnAEAYEA4AwDAgHAGAIAB4QwAAAPCGQAABoQzAAAMCGcA\nABgQzgAAMCCcAQBgQDgDAMDAKJyr6uKququq7qmqazfYf25Vfayqvl5V16zbd3pV/X5V7auqz1bV\nj23X5AEAYKfs3mxAVZ2Q5PokFyX5QpJPVNUHu/uuNcO+mOTNSS7f4BS/k+TD3f0vqmp3klMOf9oA\nALCzJnecL0hyb3ff1937k9yS5LK1A7r7ke6+I8kTa7dX1XOT/FR3v3d13BPd/ZXtmToAAOycSTif\nmeT+NesPrG6beFGSR6rqvVX1qaq6sapOPtRJAgDAoh3pNwfuTvKKJO/q7lckeTzJrx3h1wQAgG23\n6TPOSR5Mcs6a9bNWt008kOT+7v7k6voHkjzjzYVP2bt379PLS0tLWVpaGr4MAAA80/LycpaXl7fl\nXNXdBx9QtSvJ3Vl5c+BDSW5Psqe7920w9m1JHu3ut6/ZdluSX+jue1b3n9LdG/1mjt5sLgAAHF8e\neCC58MKV7zuhqtLdtZVjN73j3N0HquqqJLdm5dGOm7p7X1VdubK7b6yqM5J8MslpSZ6sqquTvKy7\nH03yliTvq6oTk3w+yRu3MlEAAFikTe847xR3nAEAvvMcS3ecfXIgAAAMCGcAABgQzgAAMCCcAQBg\nQDgDAMCAcAYAgAHhDAAAA8IZAAAGhDMAAAwIZwAAGBDOAAAwIJwBAGBAOAMAwIBwBgCAAeEMAAAD\nwhkAAAaEMwAADAhnAAAYEM4AADAgnAEAYEA4AwDAgHAGAIAB4QwAAAPCGQAABoQzAAAMCGcAABgQ\nzgAAMCCcAQBgQDgDAMCAcAYAgAHhDAAAA8IZAAAGhDMAAAwIZwAAGBDOAAAwIJwBAGBAOAMAwIBw\nBgCAAeEMAAADwhkAAAaEMwAADAhnAAAYEM4AADAgnAEAYEA4AwDAgHAGAIAB4QwAAAPCGQAABoQz\nAAAMCGcAABgQzgAAMCCcAQBgQDgDAMDAKJyr6uKququq7qmqazfYf25Vfayqvl5V12yw/4Sq+lRV\nfWg7Jg0AADtt03CuqhOSXJ/k1UnOT7Knqs5bN+yLSd6c5Def5TRXJ/mrw5gnAAAs1O7BmAuS3Nvd\n9yVJVd2S5LIkdz01oLsfSfJIVf3T9QdX1VlJLknyH5M84240AADHvgMHkv37V76++c1v/36w5Yce\nWvTM5ybhfGaS+9esP5CVmJ767SS/muT0QzgGAOA7ypNPbi08N1re6nGHc47u5DnPSU48ceXrqeWN\ntq1fvuKKRf/tz0zCecuq6rVJHu7uO6tqKUkdydcDAHjKF7+YfO1rx054Hjiw9fDcbPmp7yefvLXj\nJsu7di36n/iRNwnnB5Ocs2b9rNVtE69McmlVXZLk5CSnVdXN3f36jQbv3bv36eWlpaUsLS0NXwYA\n4Fs+8Ynkla9Mnv/87Q3P0047suFZbjFuu+Xl5SwvL2/Luaq7Dz6galeSu5NclOShJLcn2dPd+zYY\n+7Ykj3b32zfY96okv9zdlz7L6/RmcwEAmLjttuStb135DmtVVbp7S/+Jsukd5+4+UFVXJbk1K7+F\n46bu3ldVV67s7hur6owkn0xyWpInq+rqJC/r7ke3MikAADjajJ5x7u6PJDl33bYb1iw/nOTsTc5x\nWxL/3QcAwDHJJwcCAMCAcAYAgAHhDAAAA8IZAAAGhDMAAAwIZwAAGBDOAAAwIJwBAGBAOAMAwIBw\nBgCAAeEMAAADwhkAAAaEMwAADAhnAAAYEM4AADAgnAEAYEA4AwDAgHAGAIAB4QwAAAPCGQAABoQz\nAAAMCGcAABgQzgAAMCCcAQBgQDgDAMCAcAYAgAHhDAAAA8IZAAAGhDMAAAwIZwAAGBDOAAAwIJwB\nAGBAOAMAwIBwBgCAAeEMAAADwhkAAAaEMwAADAhnAAAYEM4AADAgnAEAYEA4AwDAgHAGAIAB4QwA\nAAPCGQAABoQzAAAMCGcAABgQzgAAMCCcAQBgQDgDAMCAcAYAgAHhDAAAA8IZAAAGhDMAAAyMwrmq\nLq6qu6rqnqq6doP951bVx6rq61V1zZrtZ1XVn1TVZ6vqM1X1lu2cPAAA7JTdmw2oqhOSXJ/koiRf\nSPKJqvpgd9+1ZtgXk7w5yeXrDn8iyTXdfWdVnZrkjqq6dd2xAABw1Jvccb4gyb3dfV93709yS5LL\n1g7o7ke6+46shPLa7X/b3XeuLj+aZF+SM7dl5gAAsIMm4XxmkvvXrD+QLcRvVX1vkpcn+fihHgsA\nAIu2I28OXH1M4wNJrl698wwAAMeUTZ9xTvJgknPWrJ+1um2kqnZnJZp/t7s/eLCxe/fufXp5aWkp\nS0tL05cBAIBnWF5ezvLy8racq7r74AOqdiW5OytvDnwoye1J9nT3vg3Gvi3Jo9399jXbbk7ySHdf\ns378umN7s7kAAEzcdlvy1reufIe1qirdXVs5dtM7zt19oKquSnJrVh7tuKm791XVlSu7+8aqOiPJ\nJ5OcluTJqro6ycuS/HCSn0/ymar6dJJO8uvd/ZGtTBYAABZl8qhGVkP33HXbbliz/HCSszc49E+T\n7DqcCQIAwNHAJwcCAMCAcAYAgAHhDAAAA8IZAAAGhDMAAAwIZwAAGBDOAAAwIJwBAGBAOAMAwIBw\nBgCAAeEMAAADwhkAAAaEMwAADAhnAAAYEM4AADAgnAEAYEA4AwDAgHAGAIAB4QwAAAPCGQAABoQz\nAAAMCGcAABgQzgAAMCCcAQBgQDgDAMCAcAYAgAHhDAAAA8IZAAAGhDMAAAzsXvQEAICj25NPJk88\nkRw4sPL92b4Otn+r+7Z67MMPJ7tVDtvMJQUAQ93J3/3d0R2MR+K83cmJJya7dq3E6EZfW913KMee\ncsqhHfcDP7DoK4bjjXAGgKH3vz95wxuS5z1vZ4PxpJOOzHmn+07wYCckEc4AMPb448kVVyTvec+i\nZwIsgv+GBACAAeEMAAADwhkAAAaEMwAADAhnAAAYEM4AADAgnAEAYEA4AwDAgHAGAIAB4QwAAAPC\nGQAABoQzAAAMCGcAABgQzgAAMCCcAQBgQDgDAMCAcAYAgAHhDAAAA8IZAAAGhDMAAAyMwrmqLq6q\nu6rqnqq6doP951bVx6rq61V1zaEcCwAAx4JNw7mqTkhyfZJXJzk/yZ6qOm/dsC8meXOS39zCsQAA\ncNSb3HG+IMm93X1fd+9PckuSy9YO6O5HuvuOJE8c6rEAAHAsmITzmUnuX7P+wOq2icM5FgAAjhre\nHAgAAAO7B2MeTHLOmvWzVrdNHNKxe/fufXp5aWkpS0tLw5cBAIBnWl5ezvLy8racq7r74AOqdiW5\nO8lFSR5KcnuSPd29b4Oxb0vyaHe/fQvH9mZzAYBFes97ko9+dOU7cGyqqnR3beXYTe84d/eBqroq\nya1ZebTjpu7eV1VXruzuG6vqjCSfTHJakier6uokL+vuRzc6disTBQCARZo8qpHu/kiSc9dtu2HN\n8sNJzp4eCwAAxxpvDgQAgAHhDAAAA8IZAAAGhDMAAAwIZwAAGBDOAAAwIJwBAGBAOAMAwIBwBgCA\nAeEMAAADwhkAAAaEMwAADAhnAAAYEM4AADAgnAEAYEA4AwDAgHAGAIAB4QwAAAPCGQAABoQzAAAM\nCGcAABgQzgAAMCCcAQBgQDgDAMCAcAYAgAHhDAAAA8IZAAAGhDMAAAwIZwAAGBDOAAAwIJwBAGBA\nOAMAwIBwBgCAAeEMAAADwhkAAAaEMwAADAhnAAAYEM4AADAgnAEAYEA4AwDAgHAGAIAB4QwAAAPC\nGQAABoQzAAAMCGcAABgQzgAAMCCcAQBgQDgDAMDA7kVPAAB22oEDyf79h/71mc8seubAIglnAEa6\nDy0yv/nNrcXpkTzXU+dJkhNPfPav5zzn2fddccVi/zkAi1Pdveg5JEmqqo+WuQDspMceS+6+++gI\ny4Od68CBZPfuzcPyUL+261yHcp5duxb9Tx1YlKpKd9dWjnXHGWDBfuu3khtuSM466/Ci8JRTjmyk\n7t6d1Jb+VQNwfBDOAAu2f3/yi7+YXHfdomcCwMGMfqtGVV1cVXdV1T1Vde2zjHlnVd1bVXdW1cvX\nbP+lqvrLqvqLqnpfVT1nuyYPAAA7ZdNwrqoTklyf5NVJzk+yp6rOWzfmNUle3N0vSXJlknevbn9B\nkjcneUV3/1BW7nC/blv/BAAAsAMmd5wvSHJvd9/X3fuT3JLksnVjLktyc5J098eTnF5VZ6zu25Xk\nu6tqd5JTknxhW2YOAAA7aBLOZya5f836A6vbDjbmwSRndvcXkrw9yd+sbvtyd//x1qcLAACLcUQ/\nObCqvicrd6NfmOQFSU6tKr8BEwCAY87kt2o8mOScNetnrW5bP+bsDcb8TJLPd/eXkqSq/iDJTyT5\nvY1eaO/evU8vLy0tZWlpaTA9AADY2PLycpaXl7flXJt+AEpV7Upyd5KLkjyU5PYke7p735oxlyR5\nU3e/tqouTPKO7r6wqi5IclOSH03yjSTvTfKJ7n7XBq/jA1CA70jXXZecdJJfRwewE47oB6B094Gq\nuirJrVl5tOOm7t5XVVeu7O4bu/vDVXVJVX0uyWNJ3rh67O1V9YEkn06yf/X7jVuZKAAALJKP3AZY\nMHecAXbO4dxxPqJvDgQAgOOFcAYAgAHhDAAAA8IZAAAGhDMAAAwIZwAAGBDOAAAwIJwBAGBAOAMA\nwIBwBgCAAeEMAAADwhkAAAaEMwAADAhnAAAYEM4AADAgnAEAYEA4AwDAgHAGAIAB4QwAAAPCGQAA\nBoQzAAAMCGcAABgQzgAAMCCcAQBgQDgDAMCAcAYAgAHhDAAAA8IZAAAGhDMAAAwIZwAAGBDOAAAw\nIJwBAGBAOAMAwIBwBgCAAeEMAAADwhkAAAaEMwAADAhnAAAYEM4AADAgnAEAYEA4AwDAgHAGAIAB\n4QwAAAPCGQAABoQzAAAMCGcAABgQzgAAMCCcAQBgQDgDAMCAcAYAgAHhDAAAA8IZAAAGhDMAAAyM\nwrmqLq6qu6rqnqq69lnGvLOq7q2qO6vq5Wu2n15Vv19V+6rqs1X1Y9s1eQAA2CmbhnNVnZDk+iSv\nTnJ+kj1Vdd66Ma9J8uLufkmSK5O8e83u30ny4e5+aZIfTrJvm+YOAAA7ZnLH+YIk93b3fd29P8kt\nSS5bN+ayJDcnSXd/PMnpVXVGVT03yU9193tX9z3R3V/ZvukDAMDOmITzmUnuX7P+wOq2g415cHXb\ni5I8UlXvrapPVdWNVXXy4UwYAAAW4Ui/OXB3klckeVd3vyLJ40l+7Qi/JgAAbLvdgzEPJjlnzfpZ\nq9vWjzn7Wcbc392fXF3+QJIN31yYJHv37n16eWlpKUtLS4PpAQDAxpaXl7O8vLwt56ruPviAql1J\n7k5yUZKHktyeZE9371sz5pIkb+ru11bVhUne0d0Xru67LckvdPc9VfW2JKd09zPiuap6s7kAHI+u\nuy456aSV7wAcWVWV7q6tHLvpHefuPlBVVyW5NSuPdtzU3fuq6sqV3X1jd3+4qi6pqs8leSzJG9ec\n4i1J3ldVJyb5/Lp9AABwTJg8qpHu/kiSc9dtu2Hd+lXPcuyfJ/nRrU4QAACOBj45EAAABoQzAAAM\nCGcAABgQzgAAMCCcAQBgQDgDAMCAcAYAgAHhDAAAA8IZAAAGhDMAAAwIZwAAGBDOAAAwIJwBAGBA\nOAMAwIBwBgCAAeEMAAADwhkAAAaEMwAADAhnAAAYEM4AADAgnAEAYEA4AwDAgHAGAIAB4QwAAAPC\nGQAABoQzAAAMCGcAABgQzgAAMCCcAQBgQDgDAMCAcAYAgAHhDAAAA8IZAAAGhDMAAAwIZwAAGBDO\nAAAwIJwBAGBAOAMAwIBwBgCAAeEMAAADwhkAAAaEMwAADAhnAAAYEM4AADAgnAEAYEA4AwDAgHAG\nAIAB4QwAAAPCGQAABoQzAAAMCGcAABgQzgAAMCCcAQBgYBTOVXVxVd1VVfdU1bXPMuadVXVvVd1Z\nVS9ft++EqvpUVX1oOyYNAAA7bdNwrqoTklyf5NVJzk+yp6rOWzfmNUle3N0vSXJlknevO83VSf5q\nW2YMAAALMLnjfEGSe7v7vu7en+SWJJetG3NZkpuTpLs/nuT0qjojSarqrCSXJPlv2zZrAADYYZNw\nPjPJ/WvWH1jddrAxD64Z89tJfjVJb3GOAACwcEf0zYFV9dokD3f3nUlq9QsAAI45uwdjHkxyzpr1\ns1a3rR9z9gZjfi7JpVV1SZKTk5xWVTd39+s3eqG9e/c+vby0tJSlpaXB9AAAYGPLy8tZXl7elnNV\n98GfoKiqXUnuTnJRkoeS3J5kT3fvWzPmkiRv6u7XVtWFSd7R3ReuO8+rkvxyd1/6LK/Tm80F4Hh0\n3XXJSSetfAfgyKqqdPeWnoLY9I5zdx+oqquS3JqVRztu6u59VXXlyu6+sbs/XFWXVNXnkjyW5I1b\nmQwAABytJo9qpLs/kuTcddtuWLd+1SbnuC3JbYc6QQAAOBr45EAAABgQzgAAMCCcAQBgQDgDAMCA\ncAYAgAHhDAAAA8IZAAAGRr/HGYDtsX9/8tWvfuvrK19J/vqvk5e+dNEzA2AzwhlgE9/4xrcid23w\nbmV5//7ktNNWvp773G8tv+ENi/5TArCZ6u5FzyFJUlV9tMwFOLZ1J1/72uFH7lPLybdH7uEsn3xy\nUrXYvx+A72RVle7e0k9i4QwcFbqTxx7bnru6X/1qcuKJhx+5Ty1/13ct+m8HgO1yOOHsUQ1gyw4c\nSB59dHti97HHVu7GTmL2+c8/+JjTTlsJZwDYTu44w3eYtW9OO9zg/frXk1NP3Z5HGE49Ndm1a9F/\nOwAc7zyqATzDDTck73vfM4P3qTenbccjDKeckpzgl1oCcAzxqAbwDMvLyU//dHL55d8evCed5M1p\nALAVwhmOY+edl/zIjyx6FgBwfPA/WQEAYEA4AwDAgHAGAIAB4QwAAAPCGQAABoQzAAAMCGcAABgQ\nzgAAMCCcAQBgQDgDAMCAcAYAgAHhDAAAA8IZAAAGhDMAAAwIZwAAGBDOAAAwIJwBAGBAOAMAwIBw\nBgCAAeEMAAADwhkAAAaEMwAADAhnAAAYEM4AADAgnAEAYEA4AwDAgHAGAIAB4QwAAAPCGQAABoQz\nAAAMCGcAABgQzgAAMCCcAQBgQDgDAMCAcAYAgAHhDAAAA8IZAAAGRuFcVRdX1V1VdU9VXfssY95Z\nVfdW1Z1V9fLVbWdV1Z9U1Wer6jNV9ZbtnDwAAOyUTcO5qk5Icn2SVyc5P8meqjpv3ZjXJHlxd78k\nyZVJ3r2664kk13T3+Ul+PMmb1h8LAADHgskd5wuS3Nvd93X3/iS3JLls3ZjLktycJN398SSnV9UZ\n3f233X3n6vZHk+xLcua2zR4AAHbIJJzPTHL/mvUH8sz4XT/mwfVjqup7k7w8yccPdZIAALBoO/Lm\nwKo6NckHkly9eucZAACOKbsHYx5Mcs6a9bNWt60fc/ZGY6pqd1ai+Xe7+4MHe6G9e/c+vby0tJSl\npaXB9AAAYGPLy8tZXl7elnNVdx98QNWuJHcnuSjJQ0luT7Knu/etGXNJkjd192ur6sIk7+juC1f3\n3Zzkke6+ZpPX6c3mAszt2ZNceunKdwBgRVWlu2srx256x7m7D1TVVUluzcqjHTd1976qunJld9/Y\n3R+uqkuq6nNJHkvyhtWJvTLJzyf5TFV9Okkn+fXu/shWJgsAAIsyeVQjq6F77rptN6xbv2qD4/40\nya7DmSAAABwNfHIgAAAMCGcAABgQzgAAMCCcAQBgQDgDAMCAcAYAgAHhDAAAA8IZAAAGhDMAAAwI\nZwAAGBDOAAAwIJwBAGBAOAMAwIBwBgCAAeEMAAADwhkAAAaEMwAADAhnAAAYEM4AADAgnAEAYEA4\nAwDAgHAGAIAB4QwAAAPCGQAABoQzAAAMCGcAABgQzgAAMCCcAQBgQDgDAMCAcAYAgAHhDAAAA8IZ\nAAAGhDMAAAwIZwAAGBDOAAAwIJwBAGBAOAMAwIBwBgCAAeEMAAADwhkAAAaEMwAADAhnAAAYEM4A\nADAgnAEAYEA4AwDAgHAGAIAB4QwAAAPCGQAABoQzAAAMCGcAABgQzgAAMCCcAQBgQDgDAMDAKJyr\n6uKququq7qmqa59lzDur6t6qurOqXn4oxwIAwNFu03CuqhOSXJ/k1UnOT7Knqs5bN+Y1SV7c3S9J\ncmWSd0+PhWezvLy86ClwFHJdsBHXBRtxXbDdJnecL0hyb3ff1937k9yS5LJ1Yy5LcnOSdPfHk5xe\nVWcMj4UN+YHHRlwXbMR1wUZcF2y3STifmeT+NesPrG6bjJkcCwAAR73dR+i8tZWDfvZnt3saHMvu\nvju5445Fz+LYdccdyeWXL3oWAHD8qO4++ICqC5Ps7e6LV9d/LUl3939aM+bdSf5Pd79/df2uJK9K\n8qLNjl1zjoNPBAAAtkF3b+km7+SO8yeSfH9VvTDJQ0lel2TPujEfSvKmJO9fDe0vd/fDVfXI4NjD\n+gMAAMBO2DScu/tAVV2V5NasPBN9U3fvq6orV3b3jd394aq6pKo+l+SxJG882LFH7E8DAABHyKaP\nagAAADv8yYGH80EqHL82uy6q6tyq+lhVfb2qrlnEHNl5g+viiqr689Wvj1bVDy5inuyswXVx6eo1\n8emqur2qXrmIebKzph+2VlU/WlX7q+qf7eT8WIzBz4tXVdWXq+pTq1/XbXrOnbrjvPphKPckuSjJ\nF7Ly7PTruvuuNWNek+Sq7n5tVf1Ykt/p7gt3ZIIsxPC6+EdJXpjk8iR/393/eRFzZecMr4sLk+zr\n7n+oqouz8kZkPy+OY8Pr4pTufnx1+QeT/Pfufuki5svOmFwXa8b97yRfS/Ke7v6DnZ4rO2f48+JV\nSX65uy+dnncn7zgfzgepcPza9Lro7ke6+44kTyxigizE5Lr4s+7+h9XVP4vfEf+dYHJdPL5m9dQk\nT+7g/FiM6YetvTnJB5L8v52cHAszvS4O6ZdT7GQ4b+WDVB7cYAzHFx+Sw0YO9br4t0n+6IjOiKPB\n6Lqoqsural+S/5nk3+zQ3FicTa+LqnpBksu7+79mi581wTFn+u+RH199PPh/VdXLNjvpkfoAFIAd\nUVU/nZXf5POTi54LR4fu/sMkf1hVP5nkPyT5JwueEov3jiRrn3EVzyTJHUnO6e7HVx8X/sMkP3Cw\nA3byjvODSc5Zs37W6rb1Y87eZAzHl8l1wXee0XVRVT+U5MYkl3b33+/Q3FicQ/p50d0fTfJ9VfW8\nIz0xFmpyXfzjJLdU1V8n+bkk76qq8XOtHJM2vS66+9GnHu/q7j9KcuJmPy92Mpyf/iCVqnpOVj4M\n5UPrxnwoyeuTp9/48+XufngH58jOm1wXa7lL8J1h0+uiqs5J8j+S/Ovu/r8LmCM7b3JdvHjN8iuS\nPKe7v7Sz02SHbXpddPf3rX69KCvPOf9idx/s3zUc+yY/L85Ys3xBVn5pxkF/XuzYoxqH80EqHL8m\n18Xqhf3JJKclebKqrk7ysu5+dHEz50iaXBdJ/n2S5yX5L1VVSfZ39wWLmzVH2vC6+OdV9fok38zK\nb0/4l4ubMTtheF182yE7Pkl23PC6+Lmq+ndJ9mfl58W/2uy8PgAFAAAGdvQDUAAA4FglnAEAYEA4\nAwDAgHAGAIAB4QwAAAPCGQAABoQzAAAMCGcAABj4/2kR7W431RvZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x8f32cc0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(12,10))\n",
    "plt.plot(tau,m)"
   ]
  }
 ],
 "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}