{ "cells": [ { "cell_type": "code", "execution_count": 1, "id": "938b129f-c9ac-418d-879e-f2eba52298d4", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy.integrate import quad\n", "import pastas as ps" ] }, { "cell_type": "markdown", "id": "9a50691a-85cb-487b-87ad-a238224cdda2", "metadata": {}, "source": [ "### Gamma" ] }, { "cell_type": "code", "execution_count": 2, "id": "7b8bf54b-44ac-444c-9f30-4fad27b5a546", "metadata": {}, "outputs": [], "source": [ "A = 5\n", "n = 1.5\n", "a = 50\n", "p = [A, n, a]\n", "\n", "gamma = ps.Gamma()\n", "tmax = gamma.get_tmax(p)\n", "t = np.arange(0, tmax)\n", "\n", "step = gamma.step(p)\n", "stepnum = np.zeros(len(t))\n", "for i in range(1, len(t)):\n", " stepnum[i] = quad(gamma.impulse, 0, t[i], args=(p))[0]" ] }, { "cell_type": "code", "execution_count": 3, "id": "41d2e835-bf3a-47e7-9f6b-7d0dce2d9c6b", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEGCAYAAABvtY4XAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAtAklEQVR4nO3dd3xUVfrH8c8zk04gIQFC1YReQwlVURNcFSzACqisq4u/XbGylnVd3VVX17Kr7q5l7YiigmJFLCCiBFFpht6LECDUBAgkgZSZOb8/ZkSElCHJzL2ZPO/Xa16ZuWXuN4fk4eTcO+eKMQallFKhx2F1AKWUUoGhBV4ppUKUFnillApRWuCVUipEaYFXSqkQFWZ1gBM1adLEJCcnV2vfoqIiGjRoULuBasiOmcCeueyYCeyZy46ZwJ657JgJajfX0qVL84wxTctdaYyxzSMtLc1UV2ZmZrX3DRQ7ZjLGnrnsmMkYe+ayYyZj7JnLjpmMqd1cQJapoKbqEI1SSoUoLfBKKRWitMArpVSIstVJ1vKUlZWRk5NDcXFxpdvFxcWxfv36IKXyjx0zAcTGxlJWVkZ4eLjVUZRSAWT7Ap+Tk0PDhg1JTk5GRCrcrqCggIYNGwYxWdXsmMkYQ05ODjk5OaSkpFgdRykVQAEdohGRbBFZLSIrRCSrOu9RXFxMYmJipcVd+U9EiIuLq/IvIqVU3ReMHnyGMSavJm+gxb12aXsqVT/YfohGKaUqYjwePB4PbrcLt6sMl6sMt4ThdkTidpVxNH8/u7dtwOMuw+1243GXURLVlLKIxnhKConMW4vb7QJ3GR6PG4/HTUFcR4qjmuMsPkhcbhYYN8ZjvF+NhwMJfTga3ZzIor00yV2IGA/GeDAeDxgPOU3PpTAyidjCbFrlfQfGc/xhjGFDs4tZvauUPYc/JfnAfDAeNrQZy3UX9q319hETwPngRWQbcAgwwMvGmFfK2WY8MB4gKSkpbdq0ab9YHxcXR/v27as8ltvtxul01kbsWnNipqlTp7Js2TL+85//VLj91KlTGTJkCC1atADg1ltv5dZbb6Vz5861nmvbtm0cPny4Vt+3JgoLC4mNjbU6xinsmMvKTMbjweV2U0oYpW5D+NF9mLJjGFcxJUUFRITBEWnEjoj2lLqhT/4XhLuP4vCUIcaFw1PGtrC2LIgYjMvADQXPEuYpw2lchOH9+q2jHx86huHwlDLJdQ/huAgzLsJwEY6L1z0X85z7chqbwyyJvOmUjP8qu4qX3MM5Q/YxP/KOU9bfV3YdU9wX0EW2Myvy3lPW31l6Ix95zqWvbOCDyH+csv6G0juY7elHumMFkyOeOGX91aX38r2nB5c4FvF8xLOnrB9Z8g9WmPZc4czkifCJ3mXm39ye0cGvf4OTZWRkLDXGlPu/Q6ALfEtjzG4RaQbMASYYY+ZXtH3fvn1NVtYvh+rXr19Ply5dqjyWHU9onphp8uTJZGVl8dxzz1W4fXp6Ov/+97/p27f2/yc/OVdOTo5f7Ros8+bNIz093eoYp7BjrtPNZIyh+GgRBfm5FOXvp+BoMfsbdKao1EV89hdEFmyHkkIoLcRRVkieJPJu7DUUlbi44+BDpLi2Em2KiTQlRFHKd57uXFvmLYzfRf6R1vLLEdgv3P24scxbWJdG3kCiFABQapy4COMzOZcnw28kwulgSsltiIBLwr09bwknK+YcMuMvJ9Jh+EPePzGO8F88suMGsC1xMFGmhIF73gJxIo4wcDjBGcb++N7kJ6QS6TlK9Mo3ad6iNQ6nd53DEUZB426UNTqTCE8RjQ+uRBxOHM5wHM4wcDhwxaVgohNxuoqILNiBwyGIw4mIA8SBJ7Y5EtkQp+sYzuIDOBwOxOHA4XAiDgdExeMIj8LhLkFcxb51Pz2cSHgk3377HRnnnYdD8O5TAyJSYYEP6BCNMWa37+t+EZkO9AcqLPB2NXLkSHbu3ElxcTG33XYb48ePJzY2lttuu43PPvuM6OhoZsyYQVJSEp9++imPPPIIpaWlxMfHM23aNJKSko6/V0FBAampqWzatInw8HCOHDlCamoqTz75JFlZWVx99dVER0ezcOFChg0bdrzgf/HFF/z1r3/F7XbTpEkTvv76awtbRFmprKyEnQePkltYQln2EmT/WtyF+5Fj+ThL8ilxGZ6NvZ38Y6XcceTfDPEsJFrKiPbtv9XTnOGl/wXgnfBXGehch8cIR4niqERTFtaBfEcpDSLDKIw5g100whMWgwlvAGFRlMYm80DrrkRHONmZ93dyHR4ckQ3YvnM3nbv3pGNsU75JTCY63EmYZxUlUVFEREQR4XAQAVzhe3itPuX76wpce/zVx6esH/iLV70qbat5ZUPpX+l/hm0rWRcHtKxkfQOgSSXrI4DyO51hDsHpDPzHkAJW4EWkAeAwxhT4nl8InPr3zml46NO1rNt9pNx11R2i6dqyEX+/rFul27z22mskJCRw7Ngx+vXrx6hRoygqKmLgwIE8+uij3H333UycOJH77ruPwYMHs2jRIkSE5557jieeeOIXwzINGzYkPT2dzz//nJEjRzJt2jRGjRrFmDFjeP7558vtwefm5nL99dczf/58UlJSOHjw4Gl/n8r+io8Wkrt7G4f3ZbMxsid7C0pptm0G7fO+JrrsILGuQzT25HMeZXT8+k1AeDzsFa4MmwfAURPJEWnIQWcTHI2gbZNYihoMZLmnDSa6MY6YBMIaJBAe15wZZwwiNiqMhp4+FEZHEdOgEbFOJ7FAM2DG8VQDy4t6gp9L8aF58+jYJ/2k9VE1bxhVbYHswScB031XbIQBbxtjvgjg8QLm2WefZfr06QDs3LmTzZs3ExERwaWXXgpAWloac+bMAbzX7V955ZXs2bOH4uJi2rVrd8r7/eEPf+CJJ55g5MiRvP7660ycOLHS4y9atIhzzz33+HXrCQkJtfntqSA5WnSEvds3kJ+ziZXhvdl62NA653N+degdEtx5NKaANkAbYFzxi+QRx4TobPrJTgrDE9gb242c6CYcKHHy5FldSYxrQHNpz94YJ42btiQmKoYYoDnw85mstCpS2ev8gqpdASvwxpitQM/afM/KetqBGoOfN28eX331FQsXLiQmJob09HSKi4sJDw8/frmh0+nE5XIBMGHCBO68806GDx/OzJkzeeKJU0/CnH322WRnZ/PNN9/gdrvp3r17pRmMMXppYx1RVlrC7q1r2Xg0lo2HBMeOBQzZ/TJNynbTlEPHBwQeLHmY7MhOjGrgoDCiGXkxPTENW+GMb0V00zP4oN1ZNE+MJyr8klOOMW/ePC4Z8NM7JZ2yXqmf6GWSVTh8+DCNGzcmJiaGDRs2sGjRoiq3b9WqFQBvv/12hdtde+21jB07lvvvv//4soYNG1JQUHDKtoMGDeKWW25h27Ztx4dotBdvLbfHsDW3kOzsH4ld8xaRhzaRcHQbrdy7OVPcPFZ6O7M9/bkgtpAMh4Nt8YPYEn8m4U3b0ahFR95s14v4uDjgIuDPVn87KkRpga/C0KFDeemll0hNTaVTp04MHFj5mOSDDz7ImDFjaNWqFX369CEnJ6fc7a6++mruu+8+xo4de3zZuHHjuPHGG4+fZP1J06ZNeeWVV7j88svxeDw0a9bs+JCQCrzSkmJ2blzGwS0/4Nm9grj89UwtPYe3StM5U/YyN2ISexzNyY1OYU/8EMKbd+GOzkN4Ork90RFO4BarvwVVT2mBr0JkZCSzZs06ZXlhYeHx56NHj2b06NEAjBgxghEjRgC/HDYaN24c48aNO77Pd999x+jRo4mPjz++bNSoUYwaNer463nz5h1/PmzYMIYNG1Yb35KqQt7u7azblsP8Q41ZuOoQH2a2pp2U0Q4oNNHsiGxHr3Yt6dW9J52bn0VZ4ytpHRNLa6uDK3USLfAWmDBhArNmzWLmzJlWR6n3jDHs2LSSvSu/xJmzmJYFq2lp9uFx9+Qtcy9nNoxiwRk3ENc8hWadBtAqpStdnU66Wh1cKT9ogbfA//73P6sj1Gt5u7PZtGohHx7pyvdb8nim+K8McGwgj3h2NEhlR8urad4lndWp57Dgu/mkpz9qdWSlqkULvAp5bpeLTUu/Jn/FJ7TY9w3Jnp3EGSd3OV4nrX0bDjZ9iF3t29AyuQtNavipQqXsRAu8CkmFBfl8/2M+X246TMv1k/iTeYNS42RjVCqLWl9Ok9QL+a77IO9H2JUKUVrgVcg4VlTAunnv4lj7EV2LljCj7Ca+jzyHy9sOY2liTzqcNZIe8YlWx1QqaLTAqzqt1OXhu3XbaTTnLroe+ZY0KSGXxixvNpIbBl3Msz0HEBaEOT+UsiMt8Db3ySefsG7dOu65557T3jc5OZmsrCyaNKlsQqS6afuGZfzwwyL+md2BA0UlzIjawerEi4hNu4rOAy6iaZj+aCulvwU25nK5GD58OMOHD7c6ii0cLTzMmjlv0nDd23QpW0dDE8vcdh8xun8yXdsvIjxMx9OVOpEW+CpkZ2czbNgwBg8ezIIFC2jVqhUzZsz4xVS+eXl59O3bl+zsbCZPnszHH3+M2+1m9erV3HXXXZSWlvLWW28RGRnJzJkzSUhI4Mcff+SWW24hNzeXmJgYJk6cSOfOnRk3bhwJCQksX76cPn360KNHj+PzyO/bt48bb7yRrVu3AvDiiy9y1llnlTudcSjZlX+MJZ9N4vwtj9GfInY4WrGo3W20v+B6Xmjexup4StlW3Svwr586+RLdRkKXq6D0KEwdc+r6Xr+B3ldD0QF479pfrrvu8yoPuXnzZt555x0mTpzIFVdcwYcffljp9mvWrGH58uXk5eXRq1cvHn/8cZYvX84dd9zBm2++ye2338748eN56aWX6NChA4sXL+bmm29m7ty5AGzatImvvvoKp9PJ5MmTj7/vH//4R8477zymT5+O2+0+/mna8qYzTkys2ycTjcfDxqVzeXdtIW9uiqCzODmjcT+iB99El/4XcoZezqhUlepegbdASkoKvXr1ArxTA2dnZ1e6fUZGxvEpCuLi4rjssssA6NGjB6tWraKwsJAFCxYwZszP/xmVlJQcfz5mzJhy57afO3cub775JuCdwTIuLg4ofzrjulrgjcfDyrnvEr3oKTq7NtKdIfzhnMe4dlAGreJPvT2bUqpida/AV9TjLiiAiJjKe+QNEv3qsZ8sMjLy+HOn08mxY8cICwvD4/EAUFxcXOH2Dofj+GuHw4HL5cLj8RAfH8+KFSvKj9mggd/ZKprOuK7xGEPW7Kk0XvJverm3skuSWNzlXoZdciOjYuOtjqdUnaR/51ZTcnIyS5cuBeCDDz44rX0bNWpESkoK77//PuCdD2XlypVV7nf++efz4osvAt47WB05cuS0pzO2G+PxMHP1Hu7//hjLvv2cCE8xP/R6jKS/rmHAlfcQo8VdqWrTAl9Nd9111/GTnHl5eVXvcJKpU6cyadIkevbsSbdu3ZgxY0aV+zzzzDNkZmbSo0cP0tLSWLt2LUOHDsXlcpGamsr9999f5XTGdrJ2wUw2PzaAt995A2Og5ciHaPm31fQbeQth4RFWx1OqzhNjjNUZjuvbt6/Jysr6xbL169fTpUuXKvcN1B2dasKOmcCbKycnx692DYRtaxdz5LP76XlsMftIZFO/hyiLac2QjAxL8lRm3rx5pFd60+bgs2MmsGcuO2aC2s0lIkuNMX3LW6c9eBU0h4+W8c2Lf+SM9y4i5dgaFrX9I3F3r+KcS67BobckVKrW1b2TrKrO8bjdfLA0h8dnb+Ki4jCiWvyazmP/xcBEvZ+oUoFUJwq83nS6dgVzWO7HVQtwfXI7y46dQ0rrUVw94u90axkXtOMrVZ/ZvsBHRUVx4MABEhMTtcjXAmMMhw8fJioqKqDHKSk+yrIpf6Pvzjc4LA359aCu9L9kkP4bKhVEti/wrVu3Jicnh9zc3Eq3Ky4uDnjROl12zARQVFREz549A/b+G5d/S+SnNzPIs4Mf4ofS8dpnGaDDMUoFne0LfHh4OCkpKVVuN2/ePHr37h2ERP6zYybw5goPD6/19y11efjPnI1s/PZbHo8oYuW5E+k35IpaP45Syj+2L/Cqbti5eSUffPwRLx/oz1X9LiX6ogn0jLXfJaJK1Sda4FWNGI+HH2Y8T/cVD3OtRJF61TWc36u91bGUUmiBVzVQcPggG1/9A/0LvmZtZCpNrp3M+a3bWR1LKeWjBV5Vy5Y9BwibeB693LtYmHIT/X/7CE69i5JStqK/keq0zV67lz+9t5LfOIYy4qIhDDrrYqsjKaXKoQVe+c3jdrPk9T/zztYE2rVMZ9xvH6ZlfLTVsZRSFQh4gRcRJ5AF7DLGXBro46nAKDxyiC0vXsXAY4soaz6afjfcS1S43gNVKTsLRg/+NmA90CgIx1IBkLs7m8OTLqe7axuLu97L4DF3I3rLPKVsL6C/pSLSGrgEeDWQx1GBs3nrj7hfOZ8Wrl2sTX+FAVfeo8VdqToi0L+pTwN3A54AH0cFwPdb8rj8jc3Mdgxmz+Uf0TOjnBuaK6VsK2A3/BCRS4GLjTE3i0g6cFd5Y/AiMh4YD5CUlJQ2bdq0ah2vsLCQ2NjY6gcOADtmAv9y5W9ZyHNbm3IspiV3pkWRGB3YvkBdbqtgs2MmsGcuO2aC2s2VkZFR4Q0/MMYE5AH8E8gBsoG9wFFgSmX7pKWlmerKzMys9r6BYsdMxlSda8n050zZA/Hmu8cuNvlFpbbIZBU75rJjJmPsmcuOmYyp3VxAlqmgpgasW2aMudcY09oYkwxcBcw1xvw2UMdTtWPxu/+i34q/siEqlV63TiUupvYnJVNKBYdeB6+OWzTl7wzc8jTLY86iy4QPiIpuYHUkpVQNBOVyCGPMPKPXwNvapG824tw0k6Wx6XS//WMt7kqFAL3eTfHGd1t4eNYWprZ/ip5/fI/wiEirIymlaoEW+Hpu8XtP0mH2NVzauRFPXn0WYVrclQoZWuDrsR9mvMCAdY8QE9uI/47tR7hTfxyUCiX6G11PrZz7Hr2X/Y01kb3octt0IiLtd+9YpVTNaIGvh47krKPjN7ewLawtybd8TGRUjNWRlFIBoAW+ntm0r4BXNzpZ7+xEwvgZxDZqbHUkpVSA6HXw9ci+Awf43aTlHHW0oemtX5KYoD13pUKZFvh64lhRAUdeuIjrXZ1xpl1HGy3uSoU8LfD1gMftZv0Lv6GXawtF595FvlNv1KFUfaBj8PXA4kl30KdoPks63kmv86+yOo5SKki0wIe4Hz5+nkG732BxwnAGjL3P6jhKqSDSIZoQtionnzeW5iENBtHnxlf1TkxK1TP6Gx+iDhaWcNOUZSxvcA5tJ3yi88soVQ9pgQ9BbpeLHc9dypCiz3nxt31IaBBhdSSllAV0iCYELXntTgYVL6G01yWkto63Oo5SyiLagw8xy7+cwqDdb7Ak4TL6j7rd6jhKKQtpgQ8he3duoe2Cu9nsbE/P8a9YHUcpZTEt8CHC5fYw/f0pOIyH6LGTdQIxpZQW+FDxv7lbeHx/P+YP+5LW7XtYHUcpZQNa4EPAmkWzWZI5g8v7tOKSgalWx1FK2YReRVPHHT6wj6Zf3MTjkVEkXDbB6jhKKRvRHnwdt3nyTSSYfEpGvExstN6VSSn1My3wddjSma/Tt+BrspKvp0Ovc6yOo5SyGS3wdVTevl20XXI/m8I60u+3D1sdRyllQ1rg6yBjDPfO2sXz7l8TOfplwsJ1KgKl1Km0wNdBH2VtZ86GXJpfeAdndu5jdRyllE1pga9jcndn0//zC/l9i2yuOzvF6jhKKRvTAl/H7Jx6K03NQa675DycDrE6jlLKxrTA1yHLv5xCn6JvWd72Blq37251HKWUzWmBryMKDh+k1YL72epIpu/YB6yOo5SqA7TA1xFff/gyieYQrkuf0bszKaX8ErACLyJRIrJERFaKyFoReShQxwp1S7cf4o7Nqbzc7S069km3Oo5Sqo4I5Fw0JcAQY0yhiIQD34nILGPMogAeM+S4ykp57oMvadGoCdeOOM/qOEqpOiRgPXjjVeh7Ge57mEAdL1RlffAkLx25mcfPiyI2UueGU0r5T4wJXM0VESewFGgPPG+M+Us524wHxgMkJSWlTZs2rVrHKiwsJDY2tgZpa19NMxUXHOLsrJvYHNaBw2c/hDhq5//jUGyrQLFjLjtmAnvmsmMmqN1cGRkZS40xfctdaYyp8gE0BroBbQGHP/uctH88kAl0r2y7tLQ0U12ZmZnV3jdQappp8VNXmZIHGpvtG5fXSp6fhGJbBYodc9kxkzH2zGXHTMbUbi4gy1RQUyv8m19E4oBbgLFABJALRAFJIrIIeMEYk+nP/zDGmHwRmQcMBdb4s099tyHra/rnz2Rhy2sZ1LGX1XGUUnVQZYO6HwBvAucYY/JPXCEiacA1ItLWGDOpvJ1FpClQ5ivu0cCvgMdrJ3Zoc3sM38ydTRxNSf2NzhSplKqeCgu8MeaCStYtxTu2XpkWwBu+cXgH8J4x5rNqpaxn3svayT8PnkfrK27kkobxVsdRStVRfl2WISKpQPKJ2xtjPqpsH2PMKqB3TcLVR4VHDvH17On0S+7Hxb11MjGlVPVVWeBF5DUgFVgLeHyLDVBpgVfVs/rdh3jV/TrrB89DRCcTU0pVnz89+IHGmK4BT6LYl/MjvXOmkBX3K/p21z9+lFI148+F1QtFRAt8EOx4/14EaDXqMaujKKVCgD89+DfwFvm9eKcfELwfVE0NaLJ6ZsvK7+h3eLb3ssgzO1kdRykVAvwp8K8B1wCr+XkMXtUiYwyfzFvIKFrQ7coHrY6jlAoR/hT4HcaYTwKepB6btzGXZ/d0penwz7kmPtHqOEqpEOFPgd8gIm8Dn+IdogGqvkxS+cfjdvPtp5NJSejDVQP0skilVO3xp8BH4y3sF56wTC+TrCXLv3idB4oeZWj/Zwl36v1XlFK1p8oCb4y5LhhB6iNXWSnNsv7DNkcyfS/6rdVxlFIhpsIuo4jcJyIJlawfIiKXBiZW/bDskxdoY3ZzeNBfcDidVsdRSoWYynrwq4FPRaQYWMbPs0l2AHoBXwF6wXY1FR8r4szVz7IxrBM9z7/K6jhKqRBU2WRjM4AZItIBOBvv5GFHgCnAeGPMseBEDE2fzV9CH084Zen31dqNPJRS6kT+jMFvBjYHIUu9cbTUxT+XuOjS5jWmDD7b6jhKqRClXUcLzP5yJseKjnDHhV2sjqKUCmFa4IPsWFEB52bdyuuNXyftzArPYSulVI1pgQ+ylR8/RSKHaZQ+weooSqkQV2WBF5G2IvKpiOSJyH4RmSEibYMRLtQUHy2k/eZJrInsRZcBF1kdRykV4vzpwb8NvAc0B1oC7wPvBDJUqFox4xmakI+c9xeroyil6gF/CrwYY94yxrh8jyl4pypQp6G4zM3hzQtYG9GDbmddbHUcpVQ94E+BzxSRe0QkWUTOFJG7gc9FJKGyT7qqX3ovayc3HL2ZgsunWB1FKVVP+DPZ2JW+rzectPz/8PbkdTy+CmWlJbyf+QN9z2zFgE5nWh1HKVVP+PNBJ53DtoZWzJzEByX3s6L3dL2RtlIqaPy5iibGN/HYK77XHXSSMf953G6arnqR3c5W9O8/2Oo4Sql6xJ8x+NeBUuAs3+sc4JGAJQoxq+a9T7JnBwd63aRzziilgsqfitPOGPMEUAbgm2RMxxn8FLnoWfbSlJ5D/8/qKEqpesafAl8qItH4Lo0UkXaccOs+VbH9e3fSsXQd2Z2uIzwi0uo4Sql6xp+raB4EvgDaiMhUvFMH612e/PD2rqZMcTzDtMtGWh1FKVUP+XMVzZcishQYiHdo5jZjTF7Ak9Vxm/bksyLXzW3nDyImNs7qOEqpeqjKAi8iXxtjzgc+L2eZqkD+uzfxXMQhzho0w+ooSql6qrJ7skb5PqnaREQa//TJVRFJxjsnjarAvpwf6X1oNlEN4kmI1bF3pZQ1KuvB3wDcjreYL+XnK2eOAM8HNlbdtvXzp2mCh5KOI62OopSqxyq7J+szwDMiMsEY87/TfWMRaQO8iXcWSg/wiu89Q9rRwsN03fMhKxueQ4PGza2Oo5Sqx/y5THKviDQE8H2i9SMR6ePHfi7gT8aYLnhP0N4iIl1rkLVOWD3zZeIoIuZcvaGHUspa/hT4+40xBSIyGLgIeAN4saqdjDF7jDHLfM8LgPVAq5qEtTuPx/Cv7Z14PvZWOvX9ldVxlFL1nBhT+dTuIrLcGNNbRP4JrDbGvP3TMr8P4j0xOx/obow5ctK68cB4gKSkpLRp06ad7vcAQGFhIbGxsdXat7aszHXx1NISbkiNZFDLMFtkKo8dc9kxE9gzlx0zgT1z2TET1G6ujIyMpcaYvuWuNMZU+gA+A14GfgTigUhgZVX7nbB/LN6TtJdXtW1aWpqprszMzGrvW1s+f/I6c/vD/zKlLrcxxh6ZymPHXHbMZIw9c9kxkzH2zGXHTMbUbi4gy1RQU/0ZorkCmA0MNcbkAwnAn/35n0VEwoEPganGmI/82aeu2rbuBy4u/JAr2hwh3KmTiimlrOfPJ1mPAh+d8HoPsKeq/cQ78fkkYL0x5r81CVkX5M55iuYmgi6X6MlVpZQ9BLKreTZwDTBERFb4HiF5M9KD+3fR8+CXrGoyjPgmemmkUsoe/JlsrFqMMd9RT6YV3vj5swySMppfcLvVUZRS6riAFfj6osztYVZOBIWxw7mgsz8fD1BKqeDQs4E1NGfdPt4sGojjkn9bHUUppX5BC3wNrZ37Du3ihPROzayOopRSv6AFvga2r1/Knw89xCOtFuF01IvTDUqpOkQLfA3s/fo5Skw4nYbeaHUUpZQ6hRb4aio8cohuubNYFZ9BQjOdHl8pZT9a4Ktp7RevECvHaHTOTVZHUUqpcmmBrwZjDEc2L2SLsx0d+6RbHUcppcqlBb4aFm87yPUF17P6/CmIQ5tQKWVPWp2qYdqCjcRFhzO0byeroyilVIW0wJ+mvN3ZPLL5cv6esp7oCKfVcZRSqkJa4E/T5i+eJ1aO0f+s862OopRSldICfxpcZaW02/EBq6L60rp9d6vjKKVUpbTAn4bVme/RjIO4+1xndRSllKqSFvjT4Fj+Jrk0pkfGFVZHUUqpKmmB99Pu/GP86fBovun8AGHhEVbHUUqpKul88H56L2snWzyt6H9BhtVRlFLKL9qD94Pb5eLMBX/j2jMOcEZijNVxlFLKL1rg/bDm2+n82vMlI5LdVkdRSim/aYH3gydrMgdpRI/zx1odRSml/KYFvgp5u7fTo3ABm5pfRkRklNVxlFLKb1rgq7B5zsuEiYdW599gdRSllDotehVNJTwew+KcUtxR6Qzu0NPqOEopdVq0B1+JBT8e4OmCDA4MfcHqKEopddq0wFfi+2/nkBjt4KJuza2OopRSp00LfAUO7t/FHdtv4emkmUSF67TASqm6Rwt8BTZ9+QoR4qbNeb+zOopSSlWLFvhyGI+HllvfZ0N4V5K79LU6jlJKVYsW+HKsXzybMzy7ONL1N1ZHUUqpatMCX47cJe9zhBh6XDjO6ihKKVVtASvwIvKaiOwXkTWBOkYg5B8tZXzuKF7v/CrRDRpaHUcppaotkD34ycDQAL5/QExfvosSF/zq3MFWR1FKqRoJWIE3xswHDgbq/QPBeDz0nnsN9yR+S7eWcVbHUUqpGhFjTODeXCQZ+MwYU+EdqkVkPDAeICkpKW3atGnVOlZhYSGxsbHV2vcnR3LWMXzLvXzS9AYadbu4Ru9VW5kCwY657JgJ7JnLjpnAnrnsmAlqN1dGRsZSY0z5l/sZYwL2AJKBNf5un5aWZqorMzOz2vv+ZPFTV5miB5qagsMHa/xextROpkCwYy47ZjLGnrnsmMkYe+ayYyZjajcXkGUqqKl6FY3PkfwDdD/0NWsSLyS2UWOr4yilVI1pgfdZ/+UkYqSE+MHXWx1FKaVqRSAvk3wHWAh0EpEcEfl9oI5VU8YY3s5pyrSoMXTodY7VcZRSqlYEbD54Y0ydub/d6l2HmbG/GX1H3o849I8apVRo0GoGrJs9iZ7hOYzo1dLqKEopVWvqfYEvPHKQ4Tse52+J82gUFW51HKWUqjX1vsCvm/0aMVJCnJ5cVUqFmHpf4BtveIetjmQ69km3OopSStWqel3gN6/4lg7uLeR2GqsnV5VSISdgV9HUBcuWLaGxiaPzhX+wOopSStW6etttLSpx8Y/sbjzedTpxjZtYHUcppWpdvS3ws39YT1Gpi6sGplgdRSmlAqLeDtF0n/d/vNGwKX3OuMTqKEopFRD1sgf/46oFdHRtIrLt2YiI1XGUUiog6mWBz5s/kWITThc9uaqUCmH1rsAXFeTTNXcWq+MziEtMsjqOUkoFTL0r8Ku/mERDOUajc26wOopSSgVUvSrwxhge3pnKo7F/pWOfIVbHUUqpgKpXBX7xtoOs3V9Kh/Tf6CdXlVIhr15dJln4yV/4fXRjhvccanUUpZQKuHrTjd23ayvphz7kgqQiosKdVsdRSqmAqzcFftus53BgOOOiCVZHUUqpoKgXBb6k+Cjtcz5kVcwAWqZ0sTqOUkoFRb0o8KtmTqQJ+TgH3Wx1FKWUCpqQP8nq8Rje2RJGbsRQhg2+zOo4SikVNCHfg/9mUy4fHUymZNh/9NJIpVS9EvIVb8us5+jW8CiXpra0OopSSgVVSBf4Tcu+4frDz3D/GasJd4b0t6qUUqcI6ap3bM4j5BNLt+G3Wx1FKaWCLmQL/MasufQ8toT1Kb+jYVyC1XGUUiroQrbAl3z1KIdoSI9f/9nqKEopZYmQLPDLtu5je1EYG9r9nthGja2Oo5RSlgi56+A9HsM/Zm0hJ/LPzL/iPKvjKKWUZUKuB//97HcpylnDPcM6ExMZbnUcpZSyTEj14PP27qTr4rt5OrY1XXqPtzqOUkpZKqA9eBEZKiIbRWSLiNwTyGN5PG52Tx5HA3OUmJFP4XBIIA+nlFK2F7ACLyJO4HlgGNAVGCsiXQNxLOPx4Fr6BqnFWazs9hdSug0IxGGUUqpOCWQPvj+wxRiz1RhTCkwDRgTiQBvnvMqFRTNYnDCc/qP/FIhDKKVUnSPGmMC8schoYKgx5g++19cAA4wxt5603XhgPEBSUlLatGnTTv9g7lLyN2TSqMuvcDjsc7emwsJCYmNjrY5xCjvmsmMmsGcuO2YCe+ayYyao3VwZGRlLjTF9y11pjAnIAxgDvHrC62uA/1W2T1pamqmuzMzMau8bKHbMZIw9c9kxkzH2zGXHTMbYM5cdMxlTu7mALFNBTQ3kEE0O0OaE162B3QE8nlJKqRMEssD/AHQQkRQRiQCuAj4J4PGUUkqdIGDXwRtjXCJyKzAbcAKvGWPWBup4SimlfimgH3QyxswEZgbyGEoppcoXclMVKKWU8tICr5RSIUoLvFJKhSgt8EopFaIC9knW6hCRXGB7NXdvAuTVYpzaYMdMYM9cdswE9sxlx0xgz1x2zAS1m+tMY0zT8lbYqsDXhIhkmYo+rmsRO2YCe+ayYyawZy47ZgJ75rJjJgheLh2iUUqpEKUFXimlQlQoFfhXrA5QDjtmAnvmsmMmsGcuO2YCe+ayYyYIUq6QGYNXSin1S6HUg1dKKXUCLfBKKRWi6nyBD+aNvf3Iki0iq0VkhYhk+ZYliMgcEdns+9o4wBleE5H9IrLmhGUVZhCRe31tt1FELgpyrgdFZJevvVaIyMXBzCUibUQkU0TWi8haEbnNt9yy9qokk9VtFSUiS0RkpS/XQ77lVrZVRZksbasTjuUUkeUi8pnvdfDbqqI7gdSFB95piH8E2gIRwEqgq4V5soEmJy17ArjH9/we4PEAZzgX6AOsqSoD3puhrwQigRRfWzqDmOtB4K5ytg1KLqAF0Mf3vCGwyXdsy9qrkkxWt5UAsb7n4cBiYKDFbVVRJkvb6oTj3Qm8DXzmex30tqrrPfig3di7BkYAb/ievwGMDOTBjDHzgYN+ZhgBTDPGlBhjtgFb8LZpsHJVJCi5jDF7jDHLfM8LgPVAKyxsr0oyVSRYbWWMMYW+l+G+h8HatqooU0WC9vMuIq2BS4BXTzp+UNuqrhf4VsDOE17nUPkvQ6AZ4EsRWeq7mThAkjFmD3h/eYFmFuSqKIMd2u9WEVnlG8L56U/WoOcSkWSgN95eoC3a66RMYHFb+YYcVgD7gTnGGMvbqoJMYP3P1dPA3YDnhGVBb6u6XuClnGVWXvd5tjGmDzAMuEVEzrUwiz+sbr8XgXZAL2AP8B/f8qDmEpFY4EPgdmPMkco2LWdZQHKVk8nytjLGuI0xvfDeX7m/iHSvZPOg5Kogk6VtJSKXAvuNMUv93aWcZbWSq64XeFvd2NsYs9v3dT8wHe+fWftEpAWA7+t+C6JVlMHS9jPG7PP9gnqAifz8Z2nQcolION5COtUY85FvsaXtVV4mO7TVT4wx+cA8YCg2+dk6MZMN2upsYLiIZOMdNh4iIlOwoK3qeoG3zY29RaSBiDT86TlwIbDGl+d3vs1+B8ywIF5FGT4BrhKRSBFJAToAS4IV6qcfdp9f422voOUSEQEmAeuNMf89YZVl7VVRJhu0VVMRifc9jwZ+BWzA2rYqN5PVbWWMudcY09oYk4y3Js01xvwWK9oqUGeQg/UALsZ7pcGPwN8szNEW75nwlcDan7IAicDXwGbf14QA53gH75+lZXh7Br+vLAPwN1/bbQSGBTnXW8BqYJXvh7xFMHMBg/H+KbwKWOF7XGxle1WSyeq2SgWW+46/Bnigqp/vILRVRZksbauTMqbz81U0QW8rnapAKaVCVF0folFKKVUBLfBKKRWitMArpVSI0gKvlFIhSgu8UkqFKC3wKiSJSLyI3HzC65Yi8kGAjjVSRB6oYF2h72tTEfkiEMdXqiJa4FWoigeOF3hjzG5jzOgAHetu4IXKNjDG5AJ7ROTsAGVQ6hRa4FWo+hfQzjcf+JMikiy+uehFZJyIfCwin4rINhG5VUTu9M3dvUhEEnzbtRORL3yTx30rIp1PPoiIdARKjDF5vtcpIrJQRH4QkYdP2vxj4OqAftdKnUALvApV9wA/GmN6GWP+XM767sBv8M5T8ihw1BjTG1gIXOvb5hVggjEmDbiL8nvpZwPLTnj9DPCiMaYfsPekbbOAc6r5/Sh12sKsDqCURTKNd771AhE5DHzqW74aSPXN5ngW8L53ehjAe0OGk7UAck94fTYwyvf8LeDxE9btB1rWTnylqqYFXtVXJSc895zw2oP398IB5BvvVLSVOQbEnbSsovk/onzbKxUUOkSjQlUB3lveVYvxzsG+TUTGgHeWRxHpWc6m64H2J7z+Hu8MgnDqeHtHfp7ZUKmA0wKvQpIx5gDwvYisEZEnq/k2VwO/F5GfZggt73aQ84He8vM4zm14b/byA6f27DOAz6uZRanTprNJKlVDIvIM8Kkx5qsqtpsPjDDGHApOMlXfaQ9eqZp7DIipbAMRaQr8V4u7CibtwSulVIjSHrxSSoUoLfBKKRWitMArpVSI0gKvlFIhSgu8UkqFqP8H8UmBu1hZbZsAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.plot(t[1:], step, label='analytic')\n", "plt.plot(t, stepnum, '--', label='numerical')\n", "plt.xlabel('time (d)')\n", "plt.ylabel('step (m)')\n", "plt.grid()\n", "plt.legend();" ] }, { "cell_type": "markdown", "id": "94fb2160-3637-43c7-b8f1-25d508c9b485", "metadata": {}, "source": [ "### Exponential" ] }, { "cell_type": "code", "execution_count": 4, "id": "77ca80bb-c678-46c1-aab1-f0e2ec60ebdf", "metadata": {}, "outputs": [], "source": [ "A = 5\n", "a = 50\n", "p = [A, a]\n", "\n", "exponential = ps.Exponential()\n", "tmax = exponential.get_tmax(p)\n", "t = np.arange(0, tmax)\n", "\n", "step = exponential.step(p)\n", "stepnum = np.zeros(len(t))\n", "for i in range(1, len(t)):\n", " stepnum[i] = quad(exponential.impulse, 0, t[i], args=(p))[0]" ] }, { "cell_type": "code", "execution_count": 5, "id": "6172dd00-763c-4f25-b4f0-3c34a15b3c33", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEGCAYAAABvtY4XAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAr1UlEQVR4nO3deXxU5dn/8c81k52EhAQIEMCEXZZACCoCIrgCKloBFZeW1srjzw2rVu2i1afWutT2sWptxV3RWBdAXKuSSFVACPu+70gSICvZZub+/TFDDJBlCDk5M5Pr/XrNKzNn/c4BLu7c55z7iDEGpZRSocdhdwCllFLW0AKvlFIhSgu8UkqFKC3wSikVorTAK6VUiAqzO0Bt7du3N6mpqU1at6ysjDZt2jRvIAtpXmtpXmtpXuv5mzk3N7fAGNOhzpnGmIB5ZWZmmqbKzs5u8rp20LzW0rzW0rzW8zczsNTUU1O1i0YppUKUFnillApRWuCVUipEBdRJ1rpUV1ezZ88eKioqGlwuPj6e9evXt1CqU2dn3qioKLp27Up4eLgt+1dKtYyAL/B79uwhLi6O1NRURKTe5UpKSoiLi2vBZKfGrrzGGA4ePMiePXtIS0tr8f0rpVqOpQVeRHYAJYAbcBljhp3sNioqKhot7sp/IkJSUhL5+fl2R1FKWawlWvBjjTEFp7IBLe7NS4+nUq1DwHfRKKUUgPF4cBtweQye6nJc1VUYtwu3qxqP243bQFVUezzGIIW7MVVleDwu3C4Xxfs3sn5lNGXxvXF7DNEHliHVpXg8HjAG43FTGdGOosR0PAaS9i/A6T6C8XgwxrvMkahkChKH4jGQsudjnO5KMEfneyiKSeWHxGF4PIa+u7MQjxt88zAe8mL7sTvhTMS4yNj5Chi88/CwsftUpl140h0cjRJj4XjwIrIdOIz3q/zLGPNCHctMB6YDJCcnZ2ZlZR0zPz4+nl69ejW6L7fbjdPpbI7Ylpg1axbLli3jqaeeAurOO2vWLM477zw6d+4MwG233cZtt91Gv379mj3Pli1bKCoq8nv50tJSYmNjmz2HVTTvyTPG4PIYqj2Cu/oIjsoSjKsSj6sS3FUYVxU7o/tTYcKIO7yOru7diLsKMdWIx4XD4+LD2Cm4jIOM8u/oV7kap3HhwIXT48IAf466C5fHcHXVB5zpXkYYLsKMizBcHCGKa8yfcXsMj8g/OF9yceLGgYcw3Owz7Tmn6mkA3gh/lHOca47Jv8HTjXFVjwMwO+JBMhxbjpmf6+nNpKqHAfg84l76OvYcM3+BexA/rf4NAN9E3kFXObbj4RP3mdxSfScAKyJvIkHKjpn/rms0v3bdDMDmyBsIF/cx819xXczDrp8RSRUbo6YdM+8K8xfuHNv7mGn+/p0YO3Zsbn3d31a34EcaY/aJSEfgCxHZYIxZUHsBX9F/AWDYsGFmzJgxx2xg/fr1fp2MDPSTrFFRUURERNRkrCtvVlYWw4YNo0+fPgC89tprlubJyMjwe/mcnByO/7MJZK0hr9vlorSogNLCAo6UHKayrIj8Nn0pJgbnwc203/81pqoUqSpDqktxVpfxbvwv2O1JIqN4PleVvUWkp5xIKokwVURSzZiqv7LHdOQW51zuDX/nhH1mVjzPQeK5OyyXy8LmnDD/3sLLMc4oznLsYpDJxUU4LgnDJWFUSwSO6DhinUL4kba4qhKocoTjkXA8jnAqw9owrmtXwpwCh0aztiIF4wjDOJyIw0llWDx3d+uDwyGUF9zAoqoDIE7EGQbipDqyHU+kpOMUofzg/SytLsbhdCKOMHbu3kO3vkN4vWMmYQ6h6vBzrPdUIeLA4XQg4iA5oi0fteuJQ4TKwnfYbjw4HA4QB+JwMDAilq/jOuMQ4UhxNhUOQRwORJw4RBgbHkNudDwiwpGKdeBwerfvcOBwCFOdEUwNi8QBuEwe4nDicHi3feKRbJ6/w5YWeGPMPt/PPBGZDZwJLGh4rcBzxRVXsHv3bioqKpgxYwbTp08nNjaWGTNm8NFHHxEdHc3cuXNJTk5m3rx5PPLII1RVVZGUlMSsWbNITk6u2VZJSQnp6enk5uYCUFxcTHp6Ok8++SRLly7luuuuIzo6moULFzJ+/Hj+8pe/MGzYMD777DN++9vf4na7ad++PV999ZVdh0M1M7fLRUXpYXasX0qepy157lgqCnbRZeccpKIQZ2URYVXFRLqKeSXiWha6+jK4fDH/kMeIB+JrbWtq1e9Y6BnApY6FPBvxDADlJoIjEk2FRFFsCqiKaYeJSqDAk4Y7rA0mLAoTFgVhkUzrPhBikkiuiGBJ6UAkIhpneDRhkTE4I6J5pXMmkVHRrPveTV7mw0RGxxAREUVEVDROZxiras7vXFznd51b825EnfNH1rwbVOf8c2veNfZb/ZRjPhXl5JA5akytKWNoUJdGuksS+zY8v02nhufTMr0NlhV4EWkDOIwxJb73FwH/eyrbfHjeWtbtK65zXlO7aPp3acsfLhvQ4DIvv/wyiYmJlJeXc8YZZzBp0iTKysoYPnw4f/rTn7j33nuZOXMmv//97xk1ahSLFi1CRHjxxRd54oknarplAOLi4hgzZgyff/45U6dOJSsri0mTJjFlyhSee+65moJeW35+PjfddBMLFiwgLS2NQ4cOnfT3VC3LeDwUFx7i0A/bKcnbReXB3biK97M2fBDLpT/Owu386uBDxHsKSTAljBMDS+H56pt4xz2WdNnKh5HPU24iKJVYyhyxlDvjSI4N56zERFIdQ1lYNB2JTsAZ047wNgmEx7Tl950GEd02iTbOkZQ4fkVMbDzRYWFE+3I9X5NwBHDLCbmH17zrAVxY7/fb37YtHVNSm+dgKctY2YJPBmb7rtgIA94yxnxm4f4s8/e//53Zs2cDsHv3bjZv3kxERASXXnopAJmZmXzxxReA97r9q6++mv3791NVVVXntea//OUvefTRR5k6dSqvvPIKM2fObHD/ixYtYvTo0TXbSkxMbM6vp5qgstpF3vY1HN67mfL8bbgL9xJeup8ljkG86xpNZVEe3zlvOqZ1DbCEqayP60padDSF0d3Ji8rARLcnvxw69xjAVV2H8YvkHiRGj6Uq8iaio6KJBo4OFXj6MVuru5Ws1FGWFXhjzDZgcHNus6GWtlV98Dk5OXz55ZcsXLiQmJgYxowZQ0VFBeHh4TWXGzqdTlwuFwC33347d911FxMnTiQnJ4eHHnrohG2OHDmSXbt28fXXX+N2uxk4cGCDGYwxemmjDUoKC/hh6yoK926iumAbYUW72Ojpwj+qLuFA8RHWR/ycblINQLVxctCRyOaYTvRPaUunvh1YVHQXYe1SiGnfnfiOp5HUqRt3RMdwR80eJtS8y8nJYVgQnTNQwUEvk2xEUVER7dq1IyYmhg0bNrBo0aJGl09JSQEaPkk6depUpk6dygMPPFAzLS4ujpKSkhOWPfvss7n11lvZvn17TReNtuKbh/F4KPhhNwe2raR0zzryS6vJ4kK25JXybsXN9Hb8eENYHokUxJzD2T2T6J7YjWUVTxHfoTPtu/ahfadudHI6uRa4tmaNhrv+lLKaFvhGjBs3jn/+85+kp6fTt29fhg8f3uDyDz30EFOmTCElJYXhw4ezffv2Ope76qqr+OMf/8jUqVNrpk2bNo2bb7655iTrUR06dOCFF17gyiuvxOPx0LFjx5ouIeW/6uoqdm3fyKqydqzbV0zm+sc5u/QLOlBW0wWy0ZxGWcfzGNWrA+u4h4J2bUnq3p9O3XvRMTqWCdRud/ex5Xso5S8t8I2IjIzk008/PWF6aWlpzfvJkyczefJkAC6//HIuv/zyE5afNm0a06ZNq/m8cOFCJk+eTEJCQs20SZMmMWnSpJrPOTk5Ne/Hjx/P+PHjT+GbtC7GGHZt3UDe6i8wu78nvnA93V07SAEuqnyZsLBwOrdNYkPi+dDxdGK7DiC5x2D6dOrOHMfRQVabtYdRqRanBd4Gt99+Ox9//DGffRaU55wDUllJIVtXLKB0y3e86bmIr3eU8wvPk9wd/h7FJoZdUX1Y2WEKYSnpfJ45gtTkRMKc+h+mCm1a4G3wzDPP8Oijjwb0jVmBrqLazZq1a6he8gqJeQvpWbWJdPEA8H6bbmQm96L34JvZ0f4WuvUZwsAAvstZKatogVdBwXg8bF+ziLxlH/FZcXfeOnAaPT3bmRfxGlvC+/J9yk+J7XMOqUPO5amEDt67AM85w+7YStlKC7wKWEcqq9j43w+oWvcJaYf+Sw8O0QPYEn0tPxsxmrN7DKG8yzX0i9cripSqixZ4FVCKDuezbOli3tybzH+35DPf+TvaUcrG2DPY3vMiepx9Bdd37m53TKWCghZ4ZbviokNsyHmHiA1z6H9kCYNow6bIl7jurNM40OVtOgwYyNCoGLtjKhV0tMAHuA8//JB169Zx//33n/S6qampLF26lPbt21uQ7NS4PYZvthTww5fPcMWBf3CmVHOAJJZ1uoqEM69mQca5iEOfCa/UqdACH8BcLhcTJ05k4sSJdkdpNju3rGFv9os8d2Ag35Z2ZnR0B1I7TCThrKn0yTyPZIde7aJUc9EC34gdO3Ywfvx4Ro0axXfffUdKSgpz5849ZijfgoIChg0bxo4dO3j11VeZM2cObrebNWvWcPfdd1NVVcUbb7xBZGQkn3zyCYmJiWzbto377ruP/Px8YmJimDlzJv369WPatGkkJiayfPlyhg4dyqBBg1i6dCnPPvssBw4c4Oabb2bbtm0APP/884wYMaLO4YwDiau6itVfvU3EshcZULWKrkbYmnQb119+CeedPo7IMC3qSlkh+Ar8K5ecOG3AFXD6NVB1BGZNOXH+kGsh4zooOwj//umx837+caO73Lx5M2+//TYzZ87kqquu4v33329w+TVr1rB8+XIqKiro1asXjz/+OMuXL+dXv/oVr7/+OnfeeSczZsxg5syZ9O7dm8WLF3PLLbcwf/58ADZt2sSXX36J0+nk1VdfrdnuHXfcwbnnnsvs2bNxu901d9PWNZxxUlJSo9/LaofKqnjn+51cvOAnZJjd7KMjC9NupfcFN3FDyomjbCqlmlfwFXgbpKWlMWTIEMA7NPCOHTsaXH7s2LHExcURFxdHfHw8l112GQCDBg1i1apVlJaWsnjxYqZM+fE/o8rKypr3U6ZMqXNs+/nz5/P6668D3hEs4+O9g9HWNZyxnQV+75ZVrPviVW7bewGVLkNU8mUU9h9A+nlX0yU83LZcSrU2wVfg62txl5RAREzDLfI2SX612I8XGRlZ897pdFJeXk5YWJj3gb1ARUVFvcs7HI6azw6HA5fLhcfjIT4+nhUrVtQds00bv7PVN5yxHbauWUzh548xpDibJMK4acCFTDzvHPok1/Fbl1LKcnqZQhOlpqbWPHbvvffeO6l127Zty2mnnca7774LeAfGWrlyZaPrnX/++Tz/vPeZPG63m+Li4pMeztgKG7dsYvnj4+n53kX0K/6O77tcT8nNy7hn6nj6JOtwDErZRQt8E91zzz01JzkLCgoaX+E4L774Ii+99BKDBw9mwIABzJ07t9F1nn76abKzsxk0aBCZmZmsXbuWcePG4XK5SE9P54EHHmh0OOPmtG1/Abe+tYyJL66lbfkuFnWfjvuO1Zz9P8/SQW9GUsp+xpiAeWVmZprjrVu37oRpdSkuLvZruUBhd15/j+tR2dnZNe8P7NlmFv3fdWbbA31M+gPzzF8+32AKyyqaOeGpqZ03GGheawVbXmP8zwwsNfXU1ODrg1e2qThSysqsh0nf+SoZuFneaRLzrzmbpET7r9hRSp1IC7xqlDGGDbv30/frIZxl8lkaN4Yukx7nrLR+dkdTSjUgKAq80YdONyvvb3X+2flDAb/9aAvfbYmlR9t0Dp77c4aNvMzCdEqp5hLwBT4qKoqDBw+SlJSkRb4ZGGM4ePAgUVFRDS5XXVVJ7tv/S49tb7KXJ7i+fyLnXfceYU49L69UsAj4At+1a1f27NlDfn5+g8tVVFQ0WrQCiZ15o6Ki6Nq1a73zN634Bse8Oxju3kpu7Gj+fcMI1m3YqsVdqSAT8AU+PDyctLTGb2vPyckhIyOjBRI1j0DMW+1ykfvqrxm2+1WKpC3Lz36GzIu9Qzus27DV5nRKqZOlTTIFwI6CMqb8cxEHdm5gWcJFhN+xhIyLf9r4ikqpgBXwLXhlLePxsHjOczyyIppdjm6ETfoXZw7Rm5SUCgVa4Fux0qJDbHjxRoaXzOfOuEsYMP0lOsdH2x1LKdVMtMC3Ujs3LMPx7+sZ4t7P4h63cN51f8QRpn8dlAol+i+6Ffru6/8weP4NVEgkGy6axVkjJ9gdSSllAT3J2oq4PYY/f7Ken39axrfRY3D9MpuBWtyVClmWt+BFxAksBfYaYy61en+qbmXFh1k88w7ezr+MKcNP59xL39JH5SkV4lqii2YGsB5o2wL7UnU4sHcbpS9PYrRrB38bcRHnTxxkdySlVAuwtItGRLoClwAvWrkfVb8taxbDzAvo5NrH2jEzOX/iDXZHUkq1EKtb8P8H3AvoY31ssHbhZ3T/bBoVEs2ByXMYPOhsuyMppVqQnMzIgie1YZFLgQnGmFtEZAxwT1198CIyHZgOkJycnJmVldWk/ZWWlhIbG9v0wC3M6rzL81x8sGIfj0W8TNGQW4hJ6HhK29Pjay3Na61gywv+Zx47dmyuMWZYnTPrexLIqb6APwN7gB3AD8AR4M2G1qnriU7+CrYntliZN+c/c0yv38wzlz3zX3OwtLJZtqnH11qa11rBlteY5nmik2V98MaY3xhjuhpjUoFrgPnGmOut2p/yWvr+3zjnm5/xh/Y5vHXTcBLbRNgdSSllE73RKYQsff9vDFv9ECujz2DyzX8gKlL/eJVqzVqkAhhjcoCclthXa/X9B88wbNXDrIw+g74zPiQqOsbuSEopm+mdrCHgw+9W0X/lI6yNzqDvjLla3JVSgBb4oPf52h+4c95uHk/+C73v+JCo6DZ2R1JKBQjtpA1ia7/9mK8/yyG96xXcf+PF2ueulDqGVoQgtWXVQrr/50ZuCu9AwvUP0kaLu1LqONpFE4T2bttAwgfXcERiaPOLObSL1xuFlVIn0gIfZIoP5eF580rCqaZy6nt07NrT7khKqQClBT6IVLs9vPXGCyS7D7Dn4pfp3neo3ZGUUgFMO26DhDGGB+eu4e39GXSb8AmXnH2W3ZGUUgFOC3yQ+Oadp9i00sGtYydwyeh+dsdRSgUBLfBBYO1/5zBi/SPEJI4m48Jf2R1HKRUktA8+wB3YtZmUr25jp7M7fW9+HYdD7I6klAoSWuADWGVlOYWvX0eYceG85g1i4xLsjqSUCiJa4APYZ68/SV/XRjaPeJzT+gy2O45SKshoH3yAmrtiL3duzaB8yN+55uKf2R1HKRWEtMAHoL3bNvD07OUMPe00Jk+5xO44SqkgpQU+wFRXVVL21g28IsU4r1pCmFN70ZRSTaPVI8Asee0++rg2UXDWb+ia1NbuOEqpIKYFPoCs+fYThu95lSXtLiFz/DS74yilgpwW+ABRXFhA0he3s8/RiQE3/sPuOEqpEKAFPkA89fkmFrn7ceSyfxITm2B3HKVUCNCTrAEge2Mery0/TMyYp/jJUB1nRinVPLQFb7OiwoM4sq7lgvaHuPOC3nbHUUqFEC3wNlv/2gxGeZZy3/ndiQxz2h1HKRVCtMDbaOWCOQw/PI/clOvpnXGu3XGUUiFG++Bt4q4qp+PCX7NbujD4hsfsjqOUCkHagrdJ6Zp5dDZ5FF3wFJHRsXbHUUqFIG3B22BLXin3519IQc/TuXHkBLvjKKVClLbgW5jxeHh49lLEGcHEq6fbHUcpFcK0wLew7z9+iT/vu5Gb0g7TIS7S7jhKqRCmBb4FlRYfokfun6gIT2BwWie74yilQpwW+Ba07u3fkWQK8Ux4CodDT38opaxlWYEXkSgR+V5EVorIWhF52Kp9BYN929YxZN87LEkYT5+hY+yOo5RqBaxswVcC5xljBgNDgHEiMtzC/QW0VfOexUUYqVc9ancUpVQrYVmBN16lvo/hvpexan+B7Pvth7h5/wTey3yT5JQ0u+MopVoJMca6misiTiAX6AU8Z4y5r45lpgPTAZKTkzOzsrKatK/S0lJiYwPvhiGPx8P/LcxnT3Ucfz4nmkinAIGbtz6a11qa11rBlhf8zzx27NhcY8ywOmcaYxp9Ae2AAUAPwOHPOsetnwBkAwMbWi4zM9M0VXZ2dpPXtdL3H/7TlDzY0XyVM/+Y6YGatz6a11qa11rBltcY/zMDS009NbXeSzlEJB64FZgKRAD5QBSQLCKLgH8YY7L9+Z/IGFMoIjnAOGCNP+uEgqqKI3Rd9iQ/hKUwZtRou+MopVqZhvrg3wN2A+cYY/oaY0YZY4YZY7oBjwGXi8iN9a0sIh1EJMH3Phq4ANjQfNED38o5f6Ozyad09IM4nDoUsFKqZdXbgjfGXNjAvFy8fesN6Qy85uuHdwD/NsZ81KSUQehIaRE9NvyL1RFDGDz6crvjKKVaIb/uthGRdCC19vLGmA8aWscYswrIOJVwwWzBp+8wjiLyLvg9ImJ3HKVUK9RogReRl4F0YC3g8U02QIMFvjUrOlLNvWtPIyflVR47s95fhJRSylL+tOCHG2P6W54khLwyfyXFFS5+esn5dkdRSrVi/tzotFBEtMD7qbDgB36+ZCJPdF9E/y5t7Y6jlGrF/GnBv4a3yP+Ad/gBwXujarqlyYLUptmPMowjnDlmot1RlFKtnD8F/mXgBmA1P/bBqzoUF+YzYO875MaN4Yz+Z9gdRynVyvlT4HcZYz60PEkIWDf7LwyngoSL7rc7ilJK+VXgN4jIW8A8vF00QOOXSbY2ZeUVpO38NyuihzMkvdUOmqmUCiD+FPhovIX9olrT9DLJ48xaspeXKv6XlycPtDuKUkoBfhR4Y8zPWyJIMKuocvHCgu3069WHAQNb7b1dSqkAU+9lkiLyexFJbGD+eSJyqTWxgsvyOU/zbNXvmTGqo91RlFKqRkMt+NXAPBGpAJbx42iSvfE+oelLoNU/nsjtqqb7+hcoi2hL7z6n2R1HKaVqNDTY2Fxgroj0BkbiHTysGHgTmG6MKW+ZiIFt9ZdvMsT8wNLM3yAOfYa5Uipw+NMHvxnY3AJZglJ07r/YK8kMueA6u6MopdQxtMl5CjbnZtO3ej27ev+MsPBwu+MopdQx/BouWNXtX+vCaG+mcesl/8/uKEopdQJtwTfRnsNHmL2+BM9ZNxMXX+/FRkopZZtGC7yI9BCReSJSICJ5IjJXRHq0RLhAtmr2U1zp+JppI1LtjqKUUnXypwX/FvBvoBPQBXgXeNvKUIGurKSIkTuf57qEdXRJiLY7jlJK1cmfAi/GmDeMMS7f6028QxW0Wms/fYF4KSNq9B12R1FKqXr5c5I1W0TuB7LwFvargY+P3uVqjDlkYb6AYzwekje8xiZnb/oO0yc2KaUClz8F/mrfz/85bvov8Bb8VtUfv2nJ5/T17GZR+h/1xialVEDz50antJYIEiy+WLOfw2YQ6RdPszuKUko1yJ+raGJ8A4+94Pvcu7UOMlZQWsnft3Xms8x/EdNGn7eqlAps/vQxvAJUASN8n/cAj1iWKIDNn/8fotzFXD9cBxVTSgU+fwp8T2PME0A1gG+QMbE0VQByu92MXH43s+Kfp1fHWLvjKKVUo/wp8FUiEo3v0kgR6UmtR/e1FqsWzCHFHMA9WAcVU0oFB3+uonkI+AzoJiKz8A4d3Oqe8mSWvMxh2jLg/OvtjqKUUn7x5yqa/4hILjAcb9fMDGNMgeXJAsjB/TtJL/uOJV2u5exIvXNVKRUc/LmK5itjzEFjzMfGmI+MMQUi8lVLhAsUqxfMJkw8dDnv+FsBlFIqcNXbgheRKCAGaC8i7fjxxGpbvGPStArGGP60N4M3k17jxd7pdsdRSim/NdRF8z/AnXiLeS4/Fvhi4DlrYwWOFbsOszmvlF9cOdTuKEopdVIaeibr08DTInK7MeaZk92wiHQDXsc7CqUHeMG3zaBSPu/X/F/EAc4f9IHdUZRS6qT4c5nkDyISB+C7o/UDEfGnOesC7jbGnI73BO2tItL/FLK2uPKyUgbmf0JKu2jioiPsjqOUUifFnwL/gDGmRERGARcDrwHPN7aSMWa/MWaZ730JsB5IOZWwLW31/LdoK0doc+ZP7Y6ilFInTYxpeGh3EVlujMkQkT8Dq40xbx2d5vdORFKBBcBAY0zxcfOmA9MBkpOTM7Oysk72OwBQWlpKbGzz3mEau+BBUjx72TD6BcThbNZtW5HXSprXWprXWsGWF/zPPHbs2FxjzLA6ZxpjGnwBHwH/ArYCCUAksLKx9WqtH4v3JO2VjS2bmZlpmio7O7vJ69Zl/+6txvVgvFn04l3Nut2jmjuv1TSvtTSvtYItrzH+ZwaWmnpqqj9dNFcBnwPjjDGFQCLwaz/WQ0TCgfeBWcaYoDpL+fn6Qzzr/gldx/7C7ihKKdUk/tzJegT4oNbn/cD+xtYTEQFeAtYbY/56KiHtkLWunMjOv2RGjwF2R1FKqSax8pFEI4EbgPNEZIXvNcHC/TWb7ZtW0/XAfK5M72B3FKWUajJ/BhtrEmPMNwTpsMJ5OTN5PvwNCvv80u4oSinVZPpQ0eMYj4du+z9lbfRQ2icH1VWdSil1DC3wx9mYm00Xk0dVv5/YHUUppU6JFvjjHP7+bSpNOP3GTLU7ilJKnRIt8LVUuz20yV/J+tjhxCUk2R1HKaVOiRb4Wr7dUsDEij9w8IKgu6pTKaVOoAW+lo9X7iMuMpxRg3raHUUppU6ZFnif6qoKblk3ld90WU5kWPOOO6OUUnaw7Dr4YLPxuw8ZyD4O9Uy1O4pSSjULbcH7VKyaQ4mJZsCoy+2OopRSzUILPOB2VdPr0ALWtx1JVFS03XGUUqpZaIEHNn7/OQmU4Oh/md1RlFKq2WgfPDB/l2G55yKuGKV3ryqlQkerb8F7PIY3tkXzda/7aBMXb3ccpZRqNq2+wK9dt4rOJWuZMDDZ7ihKKdWsWn0XTfE3M3k34i3K03TsGaVUaGnVLXjj8dDtwJdsjB5M23Yd7Y6jlFLNqlUX+G3rculu9lPe8xK7oyilVLNr1QU+f8l7eIyQNnKK3VGUUqrZteoCH7/3azaG96N9l9PsjqKUUs2u1Rb4vJIKrii9j0UZT9gdRSmlLNFqC3zOhnwqieDMjMF2R1FKKUu02gKf8M3D3NnmS/p3bmt3FKWUskSrLPAV5WWMKvyQ4e0KERG74yillCVaZYHfuOhTYqSS6AF6eaRSKnS1ygJfvvZjjphI+g4fb3cUpZSyTKsr8MbjoXvBN2xqk0lUdBu74yillGVa3Vg0G3btZ5e7Owm9L7U7ilJKWarVFfgvtx7hqeq7WHLBBXZHUUopS7W6LprlGzYzKCWeDnGRdkdRSilLtaoCX3w4jxfyruNXcV/ZHUUppSxnWYEXkZdFJE9E1li1j5O1feFcwsRDcv9RdkdRSinLWdmCfxUYZ+H2T1r15vkUmTb0yRhtdxSllLKcZQXeGLMAOGTV9k+W8Xjofngxm2OHER4ebnccpZSynBhjrNu4SCrwkTFmYAPLTAemAyQnJ2dmZWU1aV+lpaXExsbWO78sfyeXrL2DOe1vJmGg/Tc4NZY30Ghea2leawVbXvA/89ixY3ONMcPqnGmMsewFpAJr/F0+MzPTNFV2dnaD89+Yv8Lc/dtfmz27tjV5H82psbyBRvNaS/NaK9jyGuN/ZmCpqaemtprr4L/cUcnOdhNI6ZZmdxSllGoRreIyyYryI6Ruf5sJqXYnUUqplmPlZZJvAwuBviKyR0RutGpfjdm89CsecrzMxe322RVBKaVanGVdNMaYqVZt+2QVr/uCauOk15kBddWmUkpZqlV00XTI+5atkf2IiUu0O4pSSrWYkC/w+Qf20su1leIuenOTUqp1CfkCv2XFAhxiaJd+sd1RlFKqRYV8gf+wbACjeZEe6efYHUUppVpUyBf477YepE9aGs6wVnPJv1JKASFe4H/YuYGHiv/ApR3y7Y6ilFItLqSbtXuWfc5Y50q2pnawO4pSSrW4kG7By47/cpB40voOtTuKUkq1uJAt8MbjoVtRLjtiM3A4Q/ZrKqVUvUK28u3dto6OHKK6uz69SSnVOoVsH/zqbbvZ5+lLlyEX2h1FKaVsEbIF/pODnVgc+SiLew+2O4pSStkiJLtojMfDsi37GNEzCRGxO45SStkiJAv8zs2ryXbdwOSY5XZHUUop24Rkgd+38gsixE3PAXU/plAppVqDkCzw4bu+pYB2dO4xyO4oSillm5Ar8MbjIbV0OTvbDgXtf1dKtWIhV+D3bl1NBw5T3X2k3VGUUspWIVfgc/M8PFo9lY5DLrE7ilJK2SrkCvy3+x38O/JKUnv0tTuKUkrZKuQKvGPLF4zp6sDh0P53pVTrFlIF/uD+nTxW8UeuivjW7ihKKWW7kCrwu1bMB6BdP33AtlJKhVSBr9z+HeUmgp6DRtgdRSmlbBdSBT7p4DK2RfYjIjLS7ihKKWW7kCnwVZVHSHNto6Rjpt1RlFIqIIRMgd9UHMFFVU8gmT+3O4pSSgWE0CnwhR520IX+p/e3O4pSSgWEkHngR3reHDq0TyMuSu9gVUopCJEWfFVVNddUfcC4yLV2R1FKqYAREgW+dNdKYqWc6J46wJhSSh1laYEXkXEislFEtojI/VbtJ/GQ98lNA8++2KpdKKVU0LGswIuIE3gOGA/0B6aKiDVnQHctpDIiCeK7WbJ5pZQKRla24M8EthhjthljqoAs4HJL9lSWT1H86fqAD6WUqkWMMdZsWGQyMM4Y80vf5xuAs4wxtx233HRgOkBycnJmVlZWk/ZXVlxIm7YJp5S5JZWWlhIbG2t3DL9pXmtpXmsFW17wP/PYsWNzjTF1P4DaGGPJC5gCvFjr8w3AMw2tk5mZaZoqOzu7yevaQfNaS/NaS/Naz9/MwFJTT021sotmD1C7U7wrsM/C/SmllKrFygK/BOgtImkiEgFcA3xo4f6UUkrVYtmdrMYYl4jcBnwOOIGXjTF6J5JSSrUQS4cqMMZ8Anxi5T6UUkrVLSTuZFVKKXUiLfBKKRWitMArpVSI0gKvlFIhyrI7WZtCRPKBnU1cvT1Q0IxxrKZ5raV5raV5redv5tOMMR3qmhFQBf5UiMhSU9/tugFI81pL81pL81qvOTJrF41SSoUoLfBKKRWiQqnAv2B3gJOkea2lea2lea13yplDpg9eKaXUsUKpBa+UUqoWLfBKKRWigr7At9SDvU+FiOwQkdUiskJElvqmJYrIFyKy2feznc0ZXxaRPBFZU2tavRlF5De+Y75RRFr8aef15H1IRPb6jvMKEZkQCHlFpJuIZIvIehFZKyIzfNMD+fjWlzlQj3GUiHwvIit9eR/2TQ/IY9xA3uY9vvU9CSQYXniHId4K9AAigJVAf7tz1ZFzB9D+uGlPAPf73t8PPG5zxtHAUGBNYxnxPkR9JRAJpPn+DJwBkPch4J46lrU1L9AZGOp7Hwds8mUK5ONbX+ZAPcYCxPrehwOLgeGBeowbyNusxzfYW/At92Dv5nc58Jrv/WvAFfZFAWPMAuDQcZPry3g5kGWMqTTGbAe24P2zaDH15K2PrXmNMfuNMct870uA9UAKgX1868tcH7uPsTHGlPo+hvtehgA9xg3krU+T8gZ7gU8Bdtf6vIeG/xLaxQD/EZFc30PGAZKNMfvB+48J6GhbuvrVlzGQj/ttIrLK14Vz9NfxgMkrIqlABt4WW1Ac3+MyQ4AeYxFxisgKIA/4whgT0Me4nrzQjMc32Au81DEtEK/7HGmMGQqMB24VkdF2BzpFgXrcnwd6AkOA/cBTvukBkVdEYoH3gTuNMcUNLVrHNFuObx2ZA/YYG2PcxpgheJ//fKaIDGxg8UDN26zHN9gLfFA82NsYs8/3Mw+YjfdXqwMi0hnA9zPPvoT1qi9jQB53Y8wB3z8aDzCTH3+FtT2viITjLZSzjDEf+CYH9PGtK3MgH+OjjDGFQA4wjgA/xnBs3uY+vsFe4AP+wd4i0kZE4o6+By4C1uDN+TPfYj8D5tqTsEH1ZfwQuEZEIkUkDegNfG9DvmMc/Yfs8xO8xxlszisiArwErDfG/LXWrIA9vvVlDuBj3EFEEnzvo4ELgA0E6DGuL2+zH9+WOmts4dnoCXjP8G8Ffmd3njry9cB79nslsPZoRiAJ+ArY7PuZaHPOt/H+SliNt7VwY0MZgd/5jvlGYHyA5H0DWA2s8v2D6BwIeYFReH+dXgWs8L0mBPjxrS9zoB7jdGC5L9ca4EHf9IA8xg3kbdbjq0MVKKVUiAr2LhqllFL10AKvlFIhSgu8UkqFKC3wSikVorTAK6VUiNICr0KSiCSIyC21PncRkfcs2tcVIvJgPfNKfT87iMhnVuxfqfpogVehKgGoKfDGmH3GmMkW7ete4B8NLWCMyQf2i8hIizIodQIt8CpUPQb09I2p/aSIpIpv7HgRmSYic0RknohsF5HbROQuEVkuIotEJNG3XE8R+cw3SNx/RaTf8TsRkT5ApTGmwPc5TUQWisgSEfnjcYvPAa6z9FsrVYsWeBWq7ge2GmOGGGN+Xcf8gcC1eMf6+BNwxBiTASwEfupb5gXgdmNMJnAPdbfSRwLLan1+GnjeGHMG8MNxyy4Fzmni91HqpIXZHUApm2Qb7zjnJSJSBMzzTV8NpPtGURwBvOsdlgXwPmzheJ2B/FqfRwKTfO/fAB6vNS8P6NI88ZVqnBZ41VpV1nrvqfXZg/ffhQMoNN7hXBtSDsQfN62+8T+ifMsr1SK0i0aFqhK8j5prEuMd+3y7iEwB7+iKIjK4jkXXA71qff4W76imcGJ/ex9+HB1QKctpgVchyRhzEPhWRNaIyJNN3Mx1wI0icnQk0LoeB7kAyJAf+3Fm4H2oyxJObNmPBT5uYhalTpqOJqnUKRKRp4F5xpgvG1luAXC5MeZwyyRTrZ224JU6dY8CMQ0tICIdgL9qcVctSVvwSikVorQFr5RSIUoLvFJKhSgt8EopFaK0wCulVIjSAq+UUiHq/wNCI8XNuudvWAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.plot(t[1:], step, label='analytic')\n", "plt.plot(t, stepnum, '--', label='numerical')\n", "plt.xlabel('time (d)')\n", "plt.ylabel('step (m)')\n", "plt.grid()\n", "plt.legend();" ] }, { "cell_type": "markdown", "id": "5d007aa2-7c99-4c3c-9c2a-f413a8f5a5a3", "metadata": {}, "source": [ "### Hantush" ] }, { "cell_type": "code", "execution_count": 6, "id": "35790f67-b518-4367-90c3-2d5a425da0c2", "metadata": {}, "outputs": [], "source": [ "A = 5\n", "a = 50\n", "b = 2\n", "p = [A, a, b]\n", "\n", "hantush = ps.Hantush()\n", "tmax = hantush.get_tmax(p)\n", "t = np.arange(0, tmax)\n", "\n", "step = hantush.step(p)\n", "stepnum = np.zeros(len(t))\n", "for i in range(1, len(t)):\n", " stepnum[i] = quad(hantush.impulse, 0, t[i], args=(p))[0]" ] }, { "cell_type": "code", "execution_count": 7, "id": "e63a25d6-1b30-41a2-a02d-d4d013925c30", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEGCAYAAABvtY4XAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAsfklEQVR4nO3deXhU5fn/8fc9k42QkJAEwioJq8puQBFQQVzABVBAq7hr0YpL7c+2ttVW29rWrV+ttS64oaJxRRbFlURcAGWTVWULEohAApKFrDP374+MGCA7mTmTyf26rrkyM2f75CHcefKcM88RVcUYY0zocTkdwBhjjH9YgTfGmBBlBd4YY0KUFXhjjAlRVuCNMSZEhTkdoKqkpCRNSUlp1LZFRUW0bt26aQM1AcvVMJarYYI1FwRvtlDLtXz58lxVbVftQlUNmkdaWpo2VkZGRqO39SfL1TCWq2GCNZdq8GYLtVzAMq2hptoQjTHGhCgr8MYYE6KswBtjTIgKqpOs1SkvLyc7O5uSkpJa14uLi2PDhg0BSlV/wZgrKioKEXE6hjHGz4K+wGdnZxMbG0tKSkqtRamgoIDY2NgAJqufYMulquTl5QXlVQTGmKbl1yEaEckSkTUiskpEljVmHyUlJSQmJlqPs4mICImJibjdbqejGGP8LBA9+NGqmns0O7Di3rSsPY1pGYJ+iMYYEwK8XtRThqeiHI87Eg8uKkqL0AP78Xgq8Ho9eDwevB4Ppa07UiERcCAPKcjB66nA6/Xi9ZajHqUw4Xg8rkgiCrcTmb8NVQ94veCtqByC7HgKHsKJ3reB1vmbQb2+68JB1cuS/O7sX7WDuLxVxORvARRVL6jixcXmLhfgVSU5bymxhVmgiuJFFTwSzredL0QVuuZmEndgO5UzriuoUhbWmvWdKpf33P0+bYp3HFwGyoGwtqzpeCEAfXPeJrb0BwC+6XYpKeFN3+yifpwPXkS2AvsABZ5U1aeqWWcaMA0gOTk5LT09/ZDlcXFx9OzZs85jeTyeoBx2+CnXrFmzWLFiBQ899FCN686aNYvTTz+djh07AnDTTTdx0003ceyxxzZ5ro0bN5Kfn9/k+z1ahYWFxMTEOB3jCM0xl3o9UFGMp6wELT+At7wEKsrIC2vHj+4kpKyAY/YvQz1liLcM8ZTh8pazMnII29wpJJTlcE7RW4RpOW4tI0wrcGkFL0dMYbX0oXfFt9xc/ixu9eCmAjcewvBwh3c6X+pxnK5f8qD7McKowC0/15kLSu9hpfZiijuTB8KPKAmcVXof32lXrnYv4C/hLx6xfETJI+ygHTe5Z3N7+OtHLB9Y8hT7ieGOsFe4IWzeEct7lrxABWH8New5rgj78JBlpRpOn9KZADwU/j8muT87ZHmexpJW+iQAT4b/m7Pdh448b/e245SyRwB4KfxeRrrXHbJ8g7cr48ruA2B2xJ8Z7NqEV4ULeYDrhnZs1M/Y6NGjl6vqkOqW+bsHP0JVd4pIe+BDEflGVRdVXcFX9J8CGDJkiI4aNeqQHWzYsKFeJymD7WTmT37KFRUVRURERK0Z09PTGTJkCL179wZg5syZfsslIhze1sEgMzPTcvmUlpdTkL2ewh9zKc7Po6wgD8+BfeRE9WRj9CBKi/YzZt29tA7zElZxgAhPEZHeYl5xnc+znrOJL83h08hbj9jvPeWX85xnHL1kFx9GPnrE8oyKeD6RbgxwF9ODdZRJBGUSgYdwvBJGWEQkMa3iiPYkkV/YCXWFoRJe+dUVxqD2vegWk4Jk7WRF5MV4JRxcYagrHFxuJnQ4kbNbJZNYHMUX+9uDy42IC5fbDeLmlo6n4ImMp82BRFYUnITL5cLlciMuN7hcPJR8Et7wKKKKurL+wGREXL5lldu/kHgcLnc4YcV92Fj2a1ziQlxSeQyBv36by7BhJxFWMoDsiiJcIogI4nIhIixt0xkBXKVDyfOW4BKXbx0X4S4Xq6MTEEDKR1Dk9QCCuEDETTuXi2/CowEQz2jK8A2JSuW+e4uLLS43IoCOA5cLF/A2/vkZ82uBV9Wdvq+7RWQ2cCKwqPatgs/EiRPZvn07JSUl3HrrrUybNo2YmBhuvfVW5s+fT6tWrZgzZw7JycnMmzePv//975SVlZGYmMiTTz55SFEvKChgwIABfPfdd4SHh5Ofn8+AAQN44IEHWLZsGVOnTqVVq1YsXryYcePG8eCDDzJkyBDee+89/vjHP+LxeEhKSuLjjz92sEVMg6niLcwlP3c7ecVessO6sTu/hJTV/0dYYQ4RJXlEl+cRW7GPTNdJ/Ln8SorLytkceTlJcuhf2c9VnM3DFa1JiPBwgeRQXhFNoTuaioi2eMKiaRvfncntuxDvTuaz3NsgMgZXZAwS2Rp3ZAxnte3O2LguRLs9bC0bQURkNJFR0URFRxMV2Yr7wsO57+DRfnXEt3LvwWfDgcuOWH6i72tm5h5OHvXLI5aPPPisBzCmlkbrDAyrZXkSMKCW5XHVvrt1Zybd28UAdfWWo+pYHF/7cl+hr1EAzoX5rcCLSGvApaoFvudnAX89mn3eM28d63dWP6zQ2CGa4zu14S/n9611nWeffZaEhASKi4sZOnQokyZNoqioiGHDhnHvvffyu9/9jhkzZnDnnXcycuRIlixZgojw9NNP8/DDD/Pooz/3kmJjYxk1ahTvvPMOEydOJD09nUmTJjFlyhQee+yxgwW9qj179vDLX/6SRYsWkZqayt69exv8fRo/Ky/hwO7N7N25hdzCUtZEDSF7XzFnbPgTnQ5sIMmTSyRlxAPLPYO5tvy3AHwSMY8I8ZDvjicvLIEdrXtQ2mYwl3Q6hrbR4XxW8CBR0W2IapNIdFwSMXGJTIpL5IqoKNwuITOzzRG9vkNL4qA6grdvqhYwQcifPfhkYLbvio0w4GVVfc+Px/Ob//znP8yePRuA7du3s3HjRiIiIjjvvPMASEtL48MPK8fysrOzufjii8nJyaGsrIyuXbsesb/rrruO+++/n4kTJ/Lcc88xY8aMWo+/ZMkSTj31VFJTUwFISEhoym/P1JengtLdG8nZuZ114f3ZsqeQQWvvpc/+z0jy5hKNEg0UeI/hrrJ/ERHmol9UGJ6IXnwXfSqemE644joSkdSD148ZTPvYSNrHbKBVZBgdqxxm5CEH7RXQb9GEFr8VeFXdAgxsyn3W1tP21xh8ZmYmH330EYsXLyY6OppRo0ZRUlJCeHj4wcsN3W43FRUVANx888385je/Yfz48WRmZnLXXXcdsc8RI0aQlZXFJ598gsfjoV+/frVmUFW7tDGQVPEobM0tYt/yt2i1eQFd966nNHMHkZTTVqOZXjoDEO5qHY5EDWB1bDc0IZVWSanEdezOl116kBQTics1zunvxrRgdplkHfbv30/btm2Jjo7mm2++YcmSJXWu37lzZ6D2k6RXXHEFl1xyySG/AGJjYykoKDhi3ZNPPpnp06ezdevWg0M01otvOlqUR86GL9j33WLCcpaTXLiB0yseYW95BL8N+5BJ7sVsc3UlK+EkKhKPI7rTsbzT+2RS28UQHXGu0/GNqZEV+DqMHTuWJ554ggEDBtCnTx+GDavtpA/cfffdTJkyhc6dOzNs2DA2bdpU7XpTp07lzjvv5JJLLjn43lVXXcUNN9xw8CTrT9q1a8dTTz3FhRdeiNfrpX379geHhEwDqVKyaxMr94bz5c5yYja8wrV5D9EJ6KDCZrqwsvUwpvZqR7djUjgu+WTadojjm88+ZUwQXt1jTG2swNchMjKSBQsWHPF+YWHhweeTJ09m8uTJAEyYMIEJEyYcXPZTj/yqq67iqquuOvj+Z599xuTJk4mPjz/43qRJk5g0adLB15mZmQefjxs3jnHj7M/9xijfl032svmUbsykXe5XJHpzebHsFhboMMYk9qBThxuITDmRY/oOp0fnDvRyCaOdDm1ME7AC74Cbb76ZBQsW8O677zodJTR5Pezdm8vCbeWs/HoV934/lVQgV9uwLnIgRR2HcdnA8/nX8cfTJiocmOp0YmP8wgq8A6peNmmaiKecH1Z9wL6v0um0K5NFFf25vXw6ybGRnNjpNyQcP4p+g4Zxakyk00mNCRgr8KZZy9lfzM6376Zn1st00HyitRUrIk9E+pzPvJEj6de5DSJnOh3TGEdYgTfNTkne92xcOJN/7B3Nkqz93OLaQ2HrQZT2mUi/0y5kVFJbpyMaExSswJvmwethx9LZFC1+mh75S+iPEh/Vjl+PGcP5A//j++i5MaYqK/AmqHm9yhfLltPz/cvo7MnhB03gw6TLaHfKtTw2YBAul30AzJiaWIEPcnPnzmX9+vXccccdDd42JSWFZcuWkZSU5Idk/lWydwefL/6Me79JJmtPATOiu7O6722knX0ZY9vY7QaNqQ8r8EGsoqKC8ePHM378eKejBIwW72P1MzfSe/vr9NdoYhJn8vAlaZzWbz5hbr/eYdKYkGMFvg5ZWVmMGzeOkSNH8sUXX9C5c2fmzJlzyFS+ubm5DBkyhKysLJ5//nnefvttPB4Pa9euZfr06bhcLl588UUiIyN59913SUhIYPPmzUyfPp09e/YQHR3NjBkzOPbYY7nqqqtISEhg5cqVnHDCCfTv359ly5bx3//+l127dnHDDTewZcsWAB5//HGGDx9e7XTGzU3J/j1899ZfGZaVTgTlfNH6DGLO/ANzBp1g8/AY00jNr8A/V83cH30nwnG/gLIDMGvKkcsHXQqDp0JRHrx2xaHLrn6nzkNu3LiRV155hRkzZnDRRRfx5ptv1rr+2rVrWblyJSUlJfTs2ZP77ruPlStXctttt/HCCy/w61//mmnTpvHEE0/Qq1cvli5dyo033sjChQsB+O677/joo49wu908//zzB/d7yy23cNpppzF79mw8Hs/BT9NWN51xYmJind9XMPB6ldkrd/Dxgjd5tPxFMsJG0n7CPZwyIM3paMY0e82vwDsgNTWVQYMGAZVTA2dlZdW6/ujRo4mNjSU2NpY2bdpw/vnnA9C/f39Wr15NYWEhX3zxBVOm/PzLqLS09ODzKVOmVDu3/cKFC3nhhReAyhks4+Iqb2hQ3XTGQV/gVdm4KJ3Pln7JPXvPYGCXNFaeuoiwvXsZYMXdmCbR/Ap8TT3uggKIiK69R946sV499sNFRv786Ue3201xcTFhYWF4vV4ASkpKalzf5XIdfO1yuaioqLyBcHx8PKtWrao+Zuv6n0SsaTrjYJabvZGcV26if9ESvJLKw1NuYfzgFFwuOWT+HWPM0bGzVo2UkpLC8uXLAXjjjTcatG2bNm1ITU3l9dcrbxisqnz99dd1bjdmzBgef/xxoPIOVvn5+Q2ezthJ3opylqf/jeinR9C9cCWfpN7GMb9fysS0VLvc0Rg/sALfSLfffvvBk5y5ubkN3n7WrFk888wzDBw4kL59+zJnzpw6t3nkkUfIyMigf//+pKWlsW7dOsaOHUtFRQUDBgzgrrvuqnM6Y6ds2l3ALY+/Tf8N/8c3kQPZe+WnnHbl3bSKsrlhjPEbVQ2aR1pamh5u/fr1R7xXnfz8/HqtF2jBmmvFihUBOY7X49H35rysvf/0rg68531dsDBDvR5PjetnZGQEJFdDWa6GC9ZsoZYLWKY11NTmNwZvmo3cnO/ZMfMazi75iqs7P8g1l11O+9g67lRvjGkyVuCNXyz/ZB7dMqbTRw+wtO8f+f3kaxGXjQgaE0jNosCr3XS6SVX+VecfHq+S+eK9nLblIX5wd+TARbM56Vi77NEYJwR9gY+KiiIvL4/ExEQr8k1AVcnLy8Pj8TT5vvcVlXHrq6uI21RM+8QR9Jr2IlGxNnWvMU4J+gLfpUsXsrOz2bNnT63rlZSUEBUVfOO7wZgrKiqKoqKiJt3nxk3f8fRrs1lSNIC7J1xPv6FdbEjGGIcFfYEPDw8nNTW1zvUyMzMZPHhwABI1TLDm2rZtW5Pta9mXn9P5ncv5o5Qy9ZrFDOhxTJPt2xjTeNbFMkclc8Eb9H5nMpEuL2VTZ1txNyaIWIE3jaKqzHv5vwxfMo388HZE3LCQdr1OdDqWMaaKoB+iMcHH61X+Mncdyeu/ZEebvnSZPpfw1nYy1ZhgYwXeNIjHq9zz+ue8uHI/1596BylndkfCg+sksjGmkhV4U28VHi9vPfVXbv7hGbqNeIFrxh1nl64aE8RsDN7US4XHy6yn/sVFu/6PwqRBXDtuhBV3Y4KcFXhTJ69XmfXsI1z2w/1kJwwj9YbXISzC6VjGmDr4vcCLiFtEVorIfH8fyzQ9VeX5l1/i0uy/sStuIF1ueAtszN2YZiEQPfhbgQ0BOI7xgwc/+JYH10azvOPFdLxxLkTU/25Txhhn+bXAi0gX4FzgaX8ex/jHmxlLeS5jHRNO7MVJ1/8PiYpzOpIxpgHEnzMLisgbwD+BWOB2VT2vmnWmAdMAkpOT09LT0xt1rMLCQmJiYo4irX8011wbcvZxxoY/cSC8LbnD/4bbHZjTNc21vZwSrLkgeLOFWq7Ro0cvV9Uh1S6s6U4gR/sAzgP+53s+Cphf1zbV3dGpvkLtLi3+VluuNd/n6hd/Hq5lf0nQA99lBi6UNs/2clKw5lIN3myhlota7ujkz27ZCGC8iGQB6cDpIvKSH49nmsCu/BK+fu5WTpa1FJ/1IK16neZ0JGNMI/mtwKvqH1S1i6qmAL8AFqrqZf46njl6pRUeXn76QaZ657G339W0GX6105GMMUfBroM3B90zbz1v7unMttRLSLjgAafjGGOOUkCmKlDVTCAzEMcyjfPq4o28vHQbvxp1Mt3GWs/dmFBgc9EY1mb/SJsF00lvG8XQM992Oo4xponYEE0LV1hawQcv/JNxrqUMGHpawC6HNMb4n/1vbuGeSH+b6aXP8GOnU4kedZvTcYwxTcgKfAs2+8vNjN/8Z8oj44i/9Fmwm2QbE1Lsf3QLtfuAl+fnf0z7sAO0mvwExLRzOpIxponZSdYWyONVnl5TSo6rGyU3LiM+KcnpSMYYP7AefAv04ifrGJ7/Lvec24cOVtyNCVlW4FuYTbsLCc+4mz+Hv8gFnfY6HccY40dW4FuQCo+Xp19+mamuD9nS8Tyk8wlORzLG+JEV+Bbkhc82cs3ehznQqiM5PaY6HccY42d2krWF2PljMQUfP0hv1w70gtfw7Ix0OpIxxs+sB99C3D13HV9of/KH3IL0PtvpOMaYALAefAvw4fpdfLB+F78fey5tRvVwOo4xJkCsBx/iikor+PitZ/hPm5e47qT2TscxxgSQ9eBD3IyP1zG97FkSExMIj2jldBxjTABZgQ9h2fsOoIv/S1f3Hhj/HLjtn9uYlsSGaELYE3MXcYPrbYp7nQ+ppzgdxxgTYNalC1HLsvZywqZHCQ8Xws79h9NxjDEOsAIfgrxe5a/z1xMWdSnnnnctYfHHOB3JGOMAK/AhaO7XO1mdvZ9/X3QakQO7OB3HGOMQG4MPMeUeL5+9l85rMf9mYs9wp+MYYxxkPfgQ8+qX27i2+HlSYr24ots6HccY4yDrwYeQknIP3330LMe5vifq7LshLMLpSMYYB1mBDyEvfb6RX5a/QlFCP6TfJKfjGGMcZkM0ISK/pJwdnzxHV9ceGPe43UDbGGM9+FAx8/Ms3iw+gR3D/w49z3A6jjEmCFgPPgQUlVbwzOdbGXpsdzqfdbHTcYwxQcJ68CEgffFmHiz/B78/bo/TUYwxQcR68M1cSbmHXYue4Vr3Ski0q2aMMT+zHnwz98ZXW7nC8yb57QZDjzFOxzHGBBG/FXgRiRKRL0XkaxFZJyL3+OtYLVW5x8vWjJl0kVxiz/wDiDgdyRgTRPzZgy8FTlfVgcAgYKyIDPPj8VqcuSuzubj0TQrj+iC9znI6jjEmyPitwGulQt/LcN9D/XW8lkZVefazLbzTeiKtz7nHeu/GmCOIqv9qroi4geVAT+AxVf19NetMA6YBJCcnp6WnpzfqWIWFhcTExBxFWv/wV64NeR7u+6qEq/tFcFqXhk8q1tLa62hZroYL1myhlmv06NHLVXVItQtVtc4H0BboC3QHXPXZ5rDt44EMoF9t66WlpWljZWRkNHpbf/JXrnufeF7vu/tWLT5Q2KjtW1p7HS3L1XDBmi3UcgHLtIaaWuNlkiISB0wHLgEigD1AFJAsIkuA/6lqRn1+w6jqjyKSCYwF1tZnG1OzrblFDMt+lmGR24hy292ajDHVq20M/g1gO3CKqvZR1ZGqOkRVuwL/AiaIyLU1bSwi7UQk3ve8FXAG8E3TRW+55n+cwenuVXhPnAYR0U7HMcYEqRp78Kp6Zi3LllM5tl6bjsBM3zi8C3hNVec3KqU5aH9xOe3XP0e5K4KYEdOcjmOMCWL1+iSriAwAUqqur6pv1baNqq4GBh9NOHOkOV+sYQqLKOxzIW1bJzkdxxgTxOos8CLyLDAAWAd4fW8rUGuBN01PVflo+XqGRPTm+NG3OB3HGBPk6tODH6aqx/s9ianT4s15LNrblokXvcrxyXYzbWNM7erzQafFImIFPghkfJpJalQR5/Tv6HQUY0wzUJ8e/Ewqi/wPVE4/IFR+UHWAX5OZQ+QWlnLO1n9xXbSHqLApTscxxjQD9SnwzwKXA2v4eQzeBFhGxgdMcW1kzxCblsAYUz/1KfDfq+pcvycxNfJ6lahVMymRSNqNvNrpOMaYZqI+Bf4bEXkZmEflEA1Q92WSpuks3ZDFmIpF/JByPilRcU7HMcY0E/Up8K2oLOxV56O1yyQDaO3n8zlJyuh0xo1ORzHGNCN1FnhVtTEBBxWUlPPg9z3ZN2A2v+ua5nQcY0wzUuNlkiJyp4gk1LL8dBE5zz+xzE/eW7OD0govZwwb5HQUY0wzU1sPfg0wT0RKgBX8PJtkLyrv0PQRYFMZ+lnywt/wdEw5g7ue43QUY0wzU9tkY3OAOSLSCxhB5eRh+cBLwDRVLQ5MxJbrh53bOenAJ2zoNAmxSyONMQ1UnzH4jcDGAGQxh9ny0VMMlwqSR//K6SjGmGbInzfdNkdBvV6OyXqd9eH96NjbJuU0xjScFfggtWX5R3Tx5vDjcb9wOooxppmq13zwJvBmb2+F13Mp159+udNRjDHNVJ09eBHpLiLzRCRXRHaLyBwR6R6IcC1VhcdL+voStvS+jrj4eKfjGGOaqfoM0bwMvAZ0ADoBrwOv+DNUS7f209kMP5DJBYM7OB3FGNOM1afAi6q+qKoVvsdLVE5VYPwkdunD3BbxJqP6JDsdxRjTjNWnwGeIyB0ikiIi3UTkd8A7IpJQ2yddTeOU7tpIj+LVbGh/PpHhdorEGNN49akgF/u+Xn/Y+9dQ2ZO38fgmlJ3xNCkqtB1+pdNRjDHNXH0+6JQaiCAG8HpI2PQmi2UQJ/W3uyQaY45Ofa6iifZNPPaU73Uvm2TMP4rztpNXHsn2bhcQ7raPKBhjjk59qshzQBkw3Pc6G/i73xK1YB/tjOCM0vtIOeVSp6MYY0JAfQp8D1W9HygH8E0yZjNfNbXyEt5fuYX2sVGc2D3J6TTGmBBQnwJfJiKt8F0aKSI9qHLrPtM0ile9wb+2TuayXh7cLvv9aYw5evW5iuZu4D2gq4jMonLqYLvLUxPL/3IWZRrDiBOHOh3FGBMi6nMVzQcishwYRuXQzK2qmuv3ZC1JwS7a7VnCC+EXcmW3tk6nMcaEiPpcRfOxquap6juqOl9Vc0Xk40CEaykOrHgVF17Kjrcbexhjmk6NPXgRiQKigSQRacvPJ1bbUDknjWkiJctfYbM3hZNPGuF0FGNMCKltiOZ64NdUFvPl/Fzg84HH/BurZXko6kaKovfxf53bOB3FGBNCahyiUdVHfJ9ivV1Vu6tqqu8xUFX/W9eORaSriGSIyAYRWScitzZp8hCxp6CUV7Yn0GXw2TY8Y4xpUvW5TPIHEYkF8H2i9S0ROaEe21UA/09Vj6PyBO10EbHP31elyu43f8txbOX8gTbqZYxpWvUp8HepaoGIjATOBmYCj9e1karmqOoK3/MCYAPQ+WjChpzsr+ibNZNR8bvp0yHW6TTGmBAjqrVP7S4iK1V1sIj8E1ijqi//9F69DyKSAiwC+qlq/mHLpgHTAJKTk9PS09Mb+j0AUFhYSExMTKO29afacnX55im65nzAXzo/zdje8UGTy0mWq2GCNRcEb7ZQyzV69Ojlqjqk2oWqWusDmA88CWwG4oFI4Ou6tquyfQyVJ2kvrGvdtLQ0bayMjIxGb+tPNebyeLTwHz10wZ1jdPPugoBmUm2G7eUwy9VwwZot1HIBy7SGmlqfIZqLgPeBsar6I5AA/LY+v1lEJBx4E5ilqm/VZ5sWY/tSWpfuYXWbUXRvF3y9CWNM81efT7IeAN6q8joHyKlrO6m8JOQZYIOq/vtoQoaiPTnfU+JtR+IJ452OYowJUf6cdHwEcDlwuois8j3O8ePxmpXXS9I4pexhzhrc0+koxpgQ5bebfqrqZ9i0wtUrLWT+qh0MPqYtXROinU5jjAlRdtsgB/w4/y4e3Xs95/dPdjqKMSaEWYEPNK+XsG/nslG7cM6ALk6nMcaEMCvwAabfLyamLJdvEsfQIS7K6TjGmBBmBT7A9i17nRINJ3nIRKejGGNCnBX4QPJ6ifh2HpneQZw5qLvTaYwxIc5vV9GYI6l6+Zvrelp1bM/YmEin4xhjQpz14ANoTU4hr+7vy/FDxzgdxRjTAliBDxSvh7wF/6CHezdn9+3gdBpjTAtgBT5AvNsWM3rHk1zUcQ9x0eFOxzHGtABW4ANkz9JXKdFwOp40wekoxpgWwgp8IHg9RG9+l090MKP729UzxpjAsAIfAJ5ti4ktz2V7h7OIjbLhGWNMYFiBD4Ct337NPo2hy7ALnI5ijGlB7Dr4AHjmwKm8q91Y0i/V6SjGmBbEevB+5vFU8N7aHE47rjOtItxOxzHGtCDWg/ezpHVP80TFVvb3f9vpKMaYFsZ68P7k9dBt32LyXW047dj2TqcxxrQwVuD9qGzrF8Trj+zqMo7IMBueMcYElhV4P/ph8SuUaDjd7OoZY4wDrMD7i9dD3NYFfMoghh3Xzek0xpgWyE6y+klxaRkPlk+mdXxHznTb71FjTOBZ5fGThRt/5MWyUcR1G+h0FGNMC2UF3h+8HvZ9+hS9Y0rok2BNbIxxhlUfPziw6VMu2/Nvftl1By4Rp+MYY1ooK/B+kLM4nWKNoMfwC52OYoxpwazANzWvh6Rt77HEfQKDenR2Oo0xpgWzAt/ECr/7lDjvPvalnovLZcMzxhjnWIFvYptWfcIBjaTXyMlORzHGtHB2HXwTu7/gbApi0pib0tHpKMaYFs568E1o9/5iFm/JY/Sg3ohdPWOMcZjfCryIPCsiu0Vkrb+OEWz2vHEb/wn7D+MHWO/dGOM8f/bgnwfG+nH/wcVTQZfsd2gTFU7P5Fin0xhjjP8KvKouAvb6a//BZvfqD4jTfIr7THA6ijHGACCq6r+di6QA81W1Xy3rTAOmASQnJ6elp6c36liFhYXExMQ0atum0OqrR+hb+AUfnjiTtq2jgiZXTSxXw1iuhgvWbKGWa/To0ctVdUi1C1XVbw8gBVhb3/XT0tK0sTIyMhq97VErL9X9d3fST/454YhFjuaqheVqGMvVcMGaLdRyAcu0hppql0k2gW9z8phXNo7+Q1rOKQdjTPCzAt8E5qzbz5N6IUtPG+N0FGOMOcifl0m+AiwG+ohItohc669jOUnLS9i/4k1GdY8lKSbS6TjGGHOQ33rwqnqJv/YdTLYsncu9ZffzaefuTkcxxphD2CdZj9KB5a+xT2MYdKpdHmmMCS5W4I9C+YH99Nr3CWvjRhPbOtrpOMYYcwgr8Efh28xXiKKMVkMudTqKMcYcwQr8UcjfsJBs2jNw+NlORzHGmCNYgW+kfUVlXLn3Sl4b8AzhYW6n4xhjzBGswDfSvK93UO6Bs4cNdDqKMcZUywp8Iw3OuJw/xX9I305xTkcxxphqWYFvhO/XL6V/+Rr6HtPe6SjGGFMjK/CNsPPTFyhXN73HXOV0FGOMqZEV+AYqLyulV8581kafSFJyJ6fjGGNMjazAN9Caha+SyI9I2lVORzHGmFrZbJIN9GpWK75xj+fiUZOcjmKMMbWyHnwDbMsr4tWsaPac/GfcYeFOxzHGmFpZgW+AZe+9RJp7ExcP7ep0FGOMqZMN0dRTaVkpIzbeR5823ekQd6vTcYwxpk7Wg6+nFe+9QAfycA29xukoxhhTL1bg60FVif96BjtdHTjutIucjmOMMfViBb4eVi/5iOM837LruKsRt41qGWOaB6tW9fD5shUk0p7jz/2V01GMMaberAdfh+92FXD/jn68PXIekdE2sZgxpvmwAl+Ht959j9YRwqUn2021jTHNixX4WmzZ/C23Zd3Ak8csJKF1hNNxjDGmQazA12LH3L/jQul3jo29G2OaHyvwNVj79VcM+/Ed1nWYQHynHk7HMcaYBrMCXw2vVyl+5w+USiS9Lv6H03GMMaZRrMBXY96StbQv2caW42+kdUJHp+MYY0yj2HXwh9mVX8JdH+ykX8cZvHTByU7HMcaYRrMefBVer/L2Cw/j8hRz75ShuMIjnY5kjDGNZgW+isxZ/+T63H8y4/jVpCa1djqOMcYcFSvwPovfT+eUTQ+yPuZkhlz0R6fjGGPMUbMCD3w+7zlO+GI6O8K70f36V2xCMWNMSPBrgReRsSLyrYhsEpE7/HmsxsgvKedPr31F92V/ZXtkT5Jv/oio2LZOxzLGmCbht66qiLiBx4AzgWzgKxGZq6rr/XXM+ijI38e2tUvYt+Z9bts5hrwyFz2HPspl54wiPCrGyWjGGNOk/DkWcSKwSVW3AIhIOjABaPICn7V+KamZv2L7J+DCg6C41MtDYdeRKUMZ6FnHA94HcKuHNnKAfoBHhWs69+a086bSt5PNEmmMCT2iqv7ZschkYKyqXud7fTlwkqredNh604BpAMnJyWnp6ekNPlbR3h3ErZsJrjBUBMWF4iIj6ky2RvSio2cnZ5YsAITSiAQqYrvQqssAwqNij/r7rEthYSExMcH3l4HlahjL1XDBmi3Uco0ePXq5qg6pdqGq+uUBTAGervL6cuDR2rZJS0vTxsrIyGj0tv5kuRrGcjVMsOZSDd5soZYLWKY11FR/nmTNBrpWed0F2OnH4xljjKnCnwX+K6CXiKSKSATwC2CuH49njDGmCr+dZFXVChG5CXgfcAPPquo6fx3PGGPMofz6iR5VfRd415/HMMYYUz37JKsxxoQoK/DGGBOirMAbY0yIsgJvjDEhym+fZG0MEdkDbGvk5klAbhPGaSqWq2EsV8MEay4I3myhlqubqrarbkFQFfijISLLtKaP6zrIcjWM5WqYYM0FwZutJeWyIRpjjAlRVuCNMSZEhVKBf8rpADWwXA1juRomWHNB8GZrMblCZgzeGGPMoUKpB2+MMaYKK/DGGBOimn2BD6Ybe4tIloisEZFVIrLM916CiHwoIht9XwNyV28ReVZEdovI2irv1ZhFRP7ga8NvReTsAOe6W0R2+NptlYic40CuriKSISIbRGSdiNzqe9/RNqsll6NtJiJRIvKliHzty3WP732n26umXI7/jPmO5RaRlSIy3/fav+1V051AmsODymmINwPdgQjga+B4B/NkAUmHvXc/cIfv+R3AfQHKcipwArC2rizA8b62iwRSfW3qDmCuu4Hbq1k3kLk6Aif4nscC3/mO72ib1ZLL0TYDBIjxPQ8HlgLDgqC9asrl+M+Y73i/AV4G5vte+7W9mnsP/uCNvVW1DPjpxt7BZAIw0/d8JjAxEAdV1UXA3npmmQCkq2qpqm4FNlHZtoHKVZNA5spR1RW+5wXABqAzDrdZLblqEqhcqqqFvpfhvofifHvVlKsmAfsZE5EuwLnA04cd32/t1dwLfGdge5XX2dT+w+9vCnwgIst9NxMHSFbVHKj8zwq0dyxdzVmCoR1vEpHVviGcn/5MdSSXiKQAg6ns/QVNmx2WCxxuM99wwypgN/ChqgZFe9WQC5z/GXsY+B3grfKeX9uruRd4qeY9J6/7HKGqJwDjgOkicqqDWRrC6XZ8HOgBDAJygId87wc8l4jEAG8Cv1bV/NpWreY9v2WrJpfjbaaqHlUdROX9lk8UkX61rO50LkfbS0TOA3ar6vL6blLNew3O1dwLfFDd2FtVd/q+7gZmU/kn1S4R6Qjg+7rbqXy1ZHG0HVV1l+8/pReYwc9/igY0l4iEU1lEZ6nqW763HW+z6nIFS5v5svwIZAJjCYL2qi5XELTXCGC8iGRROZR8uoi8hJ/bq7kX+KC5sbeItBaR2J+eA2cBa315rvStdiUwx4l8PjVlmQv8QkQiRSQV6AV8GahQP/2A+1xAZbsFNJeICPAMsEFV/11lkaNtVlMup9tMRNqJSLzveSvgDOAbnG+vanM53V6q+gdV7aKqKVTWqYWqehn+bi9/nS0O1AM4h8orCzYDf3IwR3cqz3p/Daz7KQuQCHwMbPR9TQhQnleo/FO0nMrewLW1ZQH+5GvDb4FxAc71IrAGWO37we7oQK6RVP4JvBpY5Xuc43Sb1ZLL0TYDBgArfcdfC/y5rp93h3M5/jNW5Xij+PkqGr+2l01VYIwxIaq5D9EYY4ypgRV4Y4wJUVbgjTEmRFmBN8aYEGUF3hhjQpQVeBOSRCReRG6s8rqTiLzhp2NNFJE/17Cs0Pe1nYi854/jG1MTK/AmVMUDBwu8qu5U1cl+OtbvgP/VtoKq7gFyRGSEnzIYcwQr8CZU/Qvo4Zv7+wERSRHfHPQicpWIvC0i80Rkq4jcJCK/8c3TvUREEnzr9RCR93yTx30qIscefhAR6Q2Uqmqu73WqiCwWka9E5G+Hrf42MNWv37UxVViBN6HqDmCzqg5S1d9Ws7wfcCmVc5LcCxxQ1cHAYuAK3zpPATerahpwO9X30kcAK6q8fgR4XFWHAj8ctu4y4JRGfj/GNFiY0wGMcUiGVs6vXiAi+4F5vvfXAAN8szcOB16vnA4GqLz5wuE6AnuqvB4BTPI9fxG4r8qy3UCnpolvTN2swJuWqrTKc2+V114q/1+4gB+1ctrZ2hQDcYe9V9P8H1G+9Y0JCBuiMaGqgMpb3DWKVs65vlVEpkDlrI4iMrCaVTcAPau8/pzK2QLhyPH23vw8i6ExfmcF3oQkVc0DPheRtSLyQCN3MxW4VkR+miG0uttBLgIGy8/jOLdSebOXrziyZz8aeKeRWYxpMJtN0pijJCKPAPNU9aM61lsETFDVfYFJZlo668Ebc/T+AUTXtoKItAP+bcXdBJL14I0xJkRZD94YY0KUFXhjjAlRVuCNMSZEWYE3xpgQZQXeGGNC1P8HH9bqFUilPBMAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.plot(t[1:], step, label='analytic')\n", "plt.plot(t, stepnum, '--', label='numerical')\n", "plt.xlabel('time (d)')\n", "plt.ylabel('step (m)')\n", "plt.grid()\n", "plt.legend();" ] }, { "cell_type": "markdown", "id": "d5d822c8-f954-4e7c-96c6-de24895d4da8", "metadata": {}, "source": [ "### Polder - not correct yet" ] }, { "cell_type": "code", "execution_count": 8, "id": "2eaa17db-8cba-4731-926e-4f8479987063", "metadata": {}, "outputs": [], "source": [ "A = np.exp(-1)\n", "a = 100\n", "b = 0.25\n", "p = [A, a, b]\n", "\n", "polder = ps.Polder()\n", "tmax = polder.get_tmax(p)\n", "t = np.arange(0, tmax)\n", "\n", "step = polder.step(p)\n", "stepnum = np.zeros(len(t))\n", "for i in range(1, len(t)):\n", " stepnum[i] = quad(polder.impulse, 0, t[i], args=(p))[0]" ] }, { "cell_type": "code", "execution_count": 9, "id": "419cf91b-9a53-440d-958d-fc085ec0a4b8", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAtNElEQVR4nO3deXxU9b3/8dcnk4SwhB0jEDRhcQFkS0Dcg0uLtgjVUrVWW2+V2qr3tt7W2l9rl9vbe6u99da2VotW0dbKrVvFpdaNQVtEAUVl34UQhIBCmEjI9v39cSZhsk2SSU5mJvN+Ph7zmHPO93vO+XwnmfOZs32POecQERFpSVq8AxARkcSmRCEiIlEpUYiISFRKFCIiEpUShYiIRJUe7wDaa/DgwS4vLy+mecvLy+ndu3fnBpQAumO71KbkoDYlh/LyctavX7/POTcklvmTLlHk5eWxYsWKmOYNBoMUFRV1bkAJoDu2S21KDmpTcggGg8yYMeODWOfXoScREYlKiUJERKLyLVGY2QNmttfMVrdQfqWZvRd+LTWziX7FIiIisfNzj2IBMDNK+TbgHOfcBOCnwHwfYxERkRj5djLbOfeameVFKV8aMboMyPUrFhERiZ352SlgOFE865wb30q9bwMnOeeubaF8HjAPICcnp2DhwoUxxRMKhejTp09M8yay7tgutSk5qE3JIRQKMWvWrJXOucKYFuCc8+0F5AGrW6kzA1gHDGrLMgsKClysFi9eHPO8iaw7tkttSg5qU3JYvHixA1a4GLflcb2PwswmAPcDFzrn9sczFmlGbS3UVh99BTIhIwtqqiH0YXh6DdRUgauBPsdC70FQdRj2rgPnAAeu1nsNHAl9joGKg1Cy6uh0nFf32AmQnQOhUgaXLoO1ZRHltXD8GZB9LBzYCdv/4U2H8HqAMZ+CPkNg/xb4YGnT8pNnQa+BXmw7ljUtP2UuZPWFknegeEXDMoApV0FGT9jxJpS8HVEWfp/2NQikw9Yl8OF7DcstDTjFG934Iuxd07A8PQtOu8EbXbsI9m1s+Lfo0RdOnecNv/84fLy9YXmvQVB4jTe86s9QtqthefZQmPwlb3jlAigvbVjePw8mzPWG35zv/Y0iDR4N4z7nDS/9DVRVAHD89m2wZDnkjIOTLvLKX/+l938RadhkGHOBN/31X9LEiGkwssj73/nnr5uW558Fx58Ohw/Am79vWj76PMgthFAprPhD0/ITL4ShE+HgLnjnj03Lx86BY06Cj7aSt+1RcG80LJ/wBRg0Cko3wOonm84/+UvQfwTsfg/WP9u0vPCr3v928QrY9GLT8ulfh54DvP/bLYublp/5TciM302AcUsUZnYc8CRwlXNuY2v1U5ZzcOSQ90oLeBtKgE0vQ8UBqCwnd+e7sOQtGHwCjJvjlS+6CSo/gZojUF3pvZ9wIUy/3hu/57Tw9MqjiWD616HoVvjkI7hjJPUbwDrn3gZnfxsOlcCvTmka68zbveV/tA3um9G0fPbd3heqdAM8fHHT8rkPefF/+B7j1/w3rGlU/qUnvPaXvA1/vb7p/P/yopcodr4Ji25sWj68wEsU216Hv32nafno87xEseVVeOU/mpaPv9RLFBtfgH/c2bR86rVAOqx/Dt5qtDFLy4CzH/eG1/4VVj3SsDyr/9FE8f5jsG5Rw/J+I44milV/hi2vNCwfcvLRRLFygfcZNGh74dFE8eZ8L1FFGll0NFG88Rs4sKNh+UmfPZooXr8TDn8EQD7AdmDiFUcTRfB27/8t0tTrvEThamHxz2jijG8eTRTB/2panvYDL1EcKWu+PKuflyjKSyH4303L++V6iaKspPnyISd5ieLj7eR9sBAa35qWO/Vooljy86bzj5rhJYo9q2HJ7U3LT77YSxS73m6+fNIXvUSxYxm89oum5adeH9dE4ds5CjN7FCgCBgN7gB8BGQDOuXvN7H7gUo7+SapdG46fFRYWum5xZ/aBnXBwp/ePW1YC5Xu9X+SnhzdwD10Me9fC4Y+9jTjAiZ+BK/7sDf9ijDdPpAmXwyXhDdRvCrwvZaAHpIdfJ1/sLb+2Fp74qjctkOkloLR0GDnD+7JXHfY2BmnpR8vS0uG46d6X8UgI1jx5dHpawPvFfOwE78t05BBs/yeYedPr3oecBH2HQUUZfPh+RHn4NXCktyGvKGPFS49TOHVawzp9h0OPPt76y/cC5pWDN9wnx9vjORLyPjdoWN57CKRneuVHDjUt7zXI2yOoLPc+g3rhOj0HQFra0QQcWWbm/eo3835t11Y1KQ8uXe79/1VXentgjedP7+EN11Q13JOpk57ZcrkZBDKOljdhXtvaVF7dTLF5f+dG5UuWLOGcc85pWl7/uUYsPy3Ni9vVNr/++vJm2m7hv3VbylvSxvLgkiWJs53oJOE7s2M+R+HnVU9XtFJ+LdDsyetuoabaO3ywbyPs3+y9LABz7vbKH78GipcfrZ+eBflnH00UQyfCwHzoOdDbQPXI9jakda56ytvIZ/bmH2+9w5nnzjz6RQe4aWXLsaWlwdwHWy7P6Annfr/l8h59YMrVUcqz4cQoV0Zn9YW8M6KWh7JHwrEtXAPRo4/3ihZfR8oze0f/9ZbZC+jVcnlGFpDVcnndBr8ldRv8uJW3slmIKHdpgab1o81v5n0PopY3TjLtLI+mo+UpKun6ekpYh/ZA8VvecXDwEkHk4YPsYd5x2jrn3ubtKfQd5h0/zurX8J/0Uz+Nvr6IjWh1xqbWv9wiIjHS1iVWzsGulbDmKe+Y9t613vRbtnmHTwqv8ZLGkBNh4Kimv2BHntP1MYuIxECJIlbv/Mk7YRrIhONOg/N/7J2My+rnlY86N57RiYh0GiWKtqq7WiLvbO/qkBNmwuzfwUmfgZ794x2diIhvlChaU33EuwLon3d5V6oMGuNN7zMEJl8Z39hERLqAEkU0pRvg8a/Cnvdh3CVw/o9gQF68oxIR6VJKFNGUrILQHrhioXdnp4hIClKiaM6RkHeV0sTLvBvQemTHOyIRkbjRE+4a2/Gm1z3Ftte9cSUJEUlxShSRQnvh/6707oQ+5uR4RyMikhB06KmOc/D0DV4fQF9+BnoPjndEIiIJQXsUddY/53X/e94PtTchIhJBiaJO8XLIOQWmzYt3JCIiCUWJos4FP4F5i1vvWVNEJMUoUTgHpeHnJilJiIg0oUSx+124eyqsfz7ekYiIJCQlivcf8x5Tedz0eEciIpKQUjtR1NZ6D6of8ynvGRIiItJEaieKPe9D6MOjT6UTEZEmUjtRbHvNex9ZFNcwREQSWWrfmT3hchiQD32HxjsSEZGEldp7FH2GwMmfjXcUIiIJLXUTxYGd8ObvvY4ARUSkRambKHYsg7/dAuX74h2JiEhC8y1RmNkDZrbXzFa3UG5m9msz22xm75nZFL9iadbeNd79E4PHdOlqRUSSjZ97FAuAmVHKLwTGhF/zgHt8jKWpPWth8AnqtkNEpBW+JQrn3GvAR1GqzAYedp5lQH8z67rLj/auhZyxXbY6EZFkFc/LY4cDOyPGi8PTdjeuaGbz8PY6yMnJIRgMxrTCUChEMBgkUH2YMw8Wsy2UxY4Yl5VI6trVnahNyUFtSg6hUKhD88czUVgz01xzFZ1z84H5AIWFha6oqCimFQaDQernPXsXI2urGZnVL6ZlJZIG7eom1KbkoDYlh44mvngmimJgRMR4LlDSZWvP7N1lqxIRSWbxvDx2EXB1+Oqn6cBB51yTw06+2PACvPgDqKnqktWJiCQzPy+PfRR4AzjRzIrN7Ktmdr2ZXR+u8jywFdgM3Ad8w69YmtjyKqxYAGmp3YOJiEhb+LaldM5d0Uq5A27wa/1Rle2CfrlgzZ0mERGRSKl5Z3ZoD2TnxDsKEZGkkKKJYi/0PibeUYiIJIXUTBQ1VdBHiUJEpC1S82zuv6/zHoMqIiKtSs09CoC01G26iEh7pN7WsnQDPHYN7F0X70hERJJC6iWKAztgzZNwpGN9n4iIpIrUSxQVB733btDHk4hIV1CiEBGRqJQoREQkqtRLFGnp0DcXMrLiHYmISFJIvURxxr/CzWviHYWISNJIvUQhIiLtknqJ4sXbvGdRiIhIm6ReFx4734L0zHhHISKSNFJvj6LqE8joFe8oRESSRuoliuoKSNcVTyIibZV6iaLqsPYoRETaIfUSxYA86H9cvKMQEUkaqXcy+yvPxjsCEZGkknp7FCIi0i4plSistgruOxfe+0u8QxERSRoplSgCNUdg10oo3xfvUEREkkZKJYq02iPeQEbP+AYiIpJEUipRBGoqvQFdHisi0ma+Jgozm2lmG8xss5nd2kx5PzN7xszeNbM1ZnaNn/Foj0JEpP18SxRmFgDuBi4ExgJXmNnYRtVuANY65yYCRcAvzcy3jpicBSB3GvQ5xq9ViIh0O37eRzEN2Oyc2wpgZguB2cDaiDoOyDYzA/oAHwHVfgX0Se8RcO1Lfi1eRKRbMuecPws2+zww0zl3bXj8KuBU59yNEXWygUXASUA2cJlz7rlmljUPmAeQk5NTsHDhwphiCoVC9OnTJ6Z5E1l3bJfalBzUpuQQCoWYNWvWSudcYSzz+7lHYc1Ma5yVPg2sAs4FRgEvmdnrzrmyBjM5Nx+YD1BYWOiKiopiCmjVU79h0ppHYO4COOakmJaRiILBILF+JolKbUoOalNyCAaDHZrfz5PZxcCIiPFcoKRRnWuAJ51nM7ANb+/CF+nVIShdB7VVfq1CRKTb8TNRLAfGmFl++AT15XiHmSLtAM4DMLMc4ERgq18BmavxBtIy/FqFiEi349uhJ+dctZndCPwdCAAPOOfWmNn14fJ7gZ8CC8zsfbxDVd91zvl22/TRRJF6fSGKiMTK1y2mc+554PlG0+6NGC4BPuVnDJHqE0VAiUJEpK1S6s7syswBMOpcyOgd71BERJJGSv20/njgZLjkW/EOQ0QkqaTUHoWIiLRfSiWKoSUvwv+Oh4qy1iuLiAiQYokivToEB3dCWiDeoYiIJI2UShTmwt1I6T4KEZE2S7FEUesN6D4KEZE2S6lEkVZbDZYGaSnVbBGRDkmpLWZ57xEwdna8wxARSSoplSj25pzj9RwrIiJt1qaD9WY2ABgGHAa2O1d3sF9ERLq7Fvcows+z/n/hDvuWAb8H/gJ8YGaPmdmMrgqys4za/AD8ekq8wxARSSrR9igeBx4GznLOHYgsMLMC4CozG+mc+4OP8XWqQM0nUFke7zBERJJKi4nCOXdBlLKVwEpfIvKRuRpdGiuSQKqqqiguLqaioiLeodTr168f69ati3cYMcvKyiI3N5eMjM67X6yt5ygmAHmR9Z1zT3ZaFF0krbZGXYyLJJDi4mKys7PJy8vDrLmnJ3e9Q4cOkZ2dHe8wYuKcY//+/RQXF5Ofn99py211q2lmDwATgDVA3UlsByRdojBXrbuyRRJIRUVFQiWJZGdmDBo0iNLS0k5dblt+Xk93zo3t1LXGyYH+4znm2D7xDkNEIihJdC4/Ps+23Efxhpl1i0RRMvwiOO+H8Q5DRLqhBQsWcOONN7Zap6SkpH782muvZe3atX6H1mFt2aN4CC9ZfAgcwXu2tXPOTfA1MhGRbmbBggWMHz+eYcOGAXD//ffHOaK2acsexQPAVcBMYBbw2fB70pnw7o/gwYviHYaIJJg5c+ZQUFDAuHHjePDBBwHo06cP3//+95k4cSLTp09nz549ADzzzDOceuqpTJ48mfPPP79+ep1Dhw6Rn59PVVUVAGVlZeTl5fHYY4+xYsUKrrzySiZNmsThw4cpKipixYoVALzwwgtMmTKFiRMnct5553Vh61vXlj2KHc65Rb5H0gXSaqvjHYKItOAnz6xhbUnnPlRs7LC+/GjWuFbrPfDAAwwcOJDDhw9TUFDAlVdeSXl5OdOnT+dnP/sZt9xyC/fddx8/+MEPOPPMM1m2bBlmxv33388dd9zBL3/5y/plZWdnU1RUxHPPPcecOXNYuHAhl156KXPnzuXuu+/mf/7nfygsLGyw/tLSUq677jpee+018vPz+eijjzr1c+iotiSK9Wb2Z+AZvENPQHJeHqv7KESkOb/+9a956qmnANi1axebNm0iMzOTz372swAUFBTw0ksvAd4lvZdddhm7d++msrKy2ctQr732Wu644w7mzJnDgw8+yH333Rd1/cuWLePss8+uX9bAgQM7s3kd1patZk+8BPGpiGlJfHmsEoVIImrLL38/BINBXn75Zd544w169erFWWedRUVFBRkZGfVXEAUCAaqrvSMSN910EzfffDMXX3wxwWCQH//4x02WecYZZ7B9+3aWLFlCTU0N48ePjxqDcy6hr/5qdavpnLumKwLpCuZqlShEpIGDBw8yYMAAevXqxfr161m+fHmr9YcPHw7AQw891GK9q6++miuuuILbbrutflp2djaHDh1qUve0007jhhtuYNu2bfWHnhJpryJap4A/MLMWIzWzc83ss9EWbmYzzWyDmW02s1tbqFNkZqvMbI2ZLWl76O2395gz4cQL/VyFiCSZmTNnUl1dzYQJE7jtttuYOnVq1Po//vGPmTt3LmeddRaDBw9usd6VV17Jxx9/zBVXXFE/7Stf+QrXX399/cnsOkOGDGH+/PlccsklTJw4kcsuu6zjDetMzrlmX8Bs4J/AK8AvgFuAHwJ/BN4H/hcYEmX+ALAFGAlkAu8CYxvV6Q+sBY4Ljx/T0vLqXgUFBS5WixcvjnneRNYd26U2JYeOtmnt2rWdE0gnKisr65TlPPbYY+5LX/pSpyyrvRp/rosXL3bACtfK9rWlV7ROAZ8GnjazMcAZwFCgDPgTMM85d7ilecOmAZudc1sBzGxhOPlE3l3yReBJ59yO8Dr3trLMDrHaKqiuhPRMP1cjIinupptu4m9/+xvPP/98vEPpFOa8X/Kdv2CzzwMznXPXhsevAk51zt0YUedXQAYwDsgG7nLOPdzMsuYB8wBycnIKFi5cGFNMU968gSO9c1kz/nsxzZ+oQqEQffp0r65J1Kbk0NE29evXj9GjR3diRB1XU1NDIBCIdxgdsnnzZg4ePFg/HgqFmDVr1krnXGGU2Vrk55nd5k7hN85K6UABcB7e1VVvmNky59zGBjM5Nx+YD1BYWOiKiopiCqj8rTSGDBlCrPMnqmAwqDYlAbWpqXXr1iVcT63J3HtsnaysLCZPnlw/HgwGO7Q8PxNFMTAiYjwXKGmmzj7nXDlQbmavAROBjfjCgaXUY8JFRDrMz63mcmCMmeWbWSZwOdD4Du+ngbPMLN3MegGnAj4+McTR/I6OiIi0pNVEYWYjzewZM9tnZnvN7GkzG9nafM65auBG4O94G/+/OOfWmNn1ZnZ9uM464AXgPeAt4H7n3OqONChqW5z2KERE2qstW80/A38BjgWGAY8Bj7Zl4c65551zJzjnRjnnfhaedq9z7t6IOr9wzo11zo13zv2q3S1oh5JhM+HkqLd+iIh0qUWLFvHzn/88pnnz8vLYt29fJ0fUVFvOUZhz7o8R438ys+idrieo4hEXM3p8UbzDEBEBoLq6mosvvpiLL7443qFE1ZY9isVmdquZ5ZnZ8WZ2C/CcmQ2Mdud2IkqvKoOKzu2dUkSS2/bt2zn55JO57rrrGDduHLNnz27SBfi+ffvIy8sDvGdKzJkzh1mzZpGfn89vf/tb7rzzTiZPnsz06dPre37dsmULM2fOpKCggLPOOov169cD3t3ZN998MzNmzOC73/1ugwce7dmzh8997nNMnDiRiRMnsnTpUqBhN+jz58/v4k+obXsUdfeSf63R9H/BOzvc6vmKRFGw8ttwqAgu6foPWkTa4MHPNJ02bg5Muw4qP4FH5jYtn/RFmHwllO+Hv1zdsOya59q02k2bNvHoo49y3333cckll/DEE09Erb969WreeecdKioqGD16NLfffjvvvPMO3/rWt3j44Yf55je/ybx587j33nsZM2YMb775Jt/4xjd49dVXAdi4cSMvv/wygUCABQsW1C/3X//1XznnnHN46qmnqKmpIRQKAQ27QZ86dSqXXnopgwYNalPbOkNbOgVs2odu0tJVTyLSVH5+PpMmTQJg0qRJbN++PWr9GTNmkJ2dTXZ2Nv369WPWLO9ZbqeccgrvvfceoVCIpUuXMnfu0cR25Ej9UxqYO3duszf1vfrqqzz8sHfPcSAQoF+/fkDDbtB37tzJpk2bEitRhC9bvRmvP6Z54S49TnTOPet7dJ3MHJDAXfmKpLxoewCZvaKX9x7U5j2Ixnr06FE/HAgEqKqqIj09ndraWgAqKiparJ+WllY/npaWRnV1NbW1tfTv359Vq1Y1H2rv3m2OrXE36EVFRU3i8VtbzlE8CFQCp4fHi4H/9C0iX+nyWBFpm7y8PFauXAnA448/3q55+/btS35+Po899hjgdb767rvvtjrfeeedxz333AN4XYmUlZU16QZ92bJl7WxJx7VlqznKOXcHUAUQ7gwwSX+W15K0oYtIl/r2t7/NPffcw+mnnx7TJaiPPPIIf/jDH5g4cSLjxo3j6aefbnWeu+66i8WLF3PKKadQUFDAmjVrmnSDPn369Fia0zGtdS8LLMXrh+nt8Pgo4K1Yu6vt6Ksj3Yxv+NN3nNvw95jnT1Tqvjo5qE1NdeduxuOps7sZb8sexY/x7p4eYWaP4D2f4rt+JC2/lQy/CE74VOsVRUSkXluuenrRzFYC0/GO2/ybc87/WwF9kHV4D5Tvg94tP5VKREQaaktfT6845/Y7555zzj3rnNtnZq90RXCdbcrbt8CrP413GCIiSaXFPQozywJ6AYPNbABHzwL3xevzKQnpPgqRROOcw3TZeqdxPjyMLtqhp68B38RLCis5uoUtA+7u9Ei6iv4hRRJGVlYW+/fvZ9CgQUoWncA5x/79+8nKyurU5UZ7ZvZdwF1mdpNz7jeduta40X0UIokkNzeX4uJiSktL4x1KvYqKik7f0HalrKwscnNzO3WZbenr6UMzy3bOHTKzHwBTgP90zr3dqZF0AXM69CSSSDIyMsjPT6xegoLBYIPHiErbbri7LZwkzgQ+DTwE3ONvWP7YMurLMO5z8Q5DRCSptCVR1ITfPwPc45x7Gsj0LyT/fDj0Asg7I95hiIgklbYkil1m9nvgC8DzZtajjfMlnN6h7VBWEu8wRESSSls2+F/Ae+71TOfcAWAg8B0/g/LL5HduhaXd5Ly8iEgXacud2Z8AT0aM7wZ2+xmUf3TVk4hIe6XUVtN8uBFFRKS7S6lEoT0KEZH2S6mtpjmnO7NFRNoppRLF+pNu1H0UIiLtlFKJYm9OEQzTHZciIu3ha6Iws5lmtsHMNpvZrVHqTTWzGjP7vJ/x9DuwFg7s8HMVIiLdjm+JwswCeL3MXgiMBa4ws7Et1Lsd714NX01a9X1Y+ZDfqxER6Vb83KOYBmx2zm11zlUCC4HZzdS7CXgC2OtjLGE6mS0i0l5t6T02VsOBnRHjxcCpkRXMbDjwOeBcYGpLCzKzecA8gJycHILBYEwBFeHY/sEOtsc4f6IKhUIxfyaJSm1KDmpTcgiFQh2a389E0dxP98Z3vP0K+K5zribaQ0ucc/OB+QCFhYWuqKio/dE4B0HIy8snL5b5E1gwGCSmzySBqU3JQW1KDh1NfH4mimJgRMR4LtC4R75CYGE4SQwGLjKzaufcXzs9mrq7snXoSUSkXfxMFMuBMWaWD+wCLge+GFnBOVf/xBIzWwA860uS8FbA6nG3Ml73UYiItItvicI5V21mN+JdzRQAHnDOrTGz68Pl9/q17maZsW/IaTDkxC5drYhIsvNzjwLn3PPA842mNZsgnHNf8TMWamsZuH8lfHQcDBzp66pERLqT1Lkzu7aKCe//B6x+svW6IiJSL3UShU5mi4jEJHUSRd2VuepmXESkXVJnq+lqwwPaoxARaY8UShQ69CQiEovUSRTpWbw74Scwdk68IxERSSqpkygC6Xw8cBIMOD7ekYiIJJXUSRTVlQzZ+w/YvyXekYiIJJXUSRSVIcat/QVsejHekYiIJJXUSRT1dDJbRKQ9UidRON1HISISi9TZatbdR6HLY0VE2iV1EkWTZyaJiEhbpE6i6DmAtyffASdfHO9IRESSSuokikAGZf1OhOyceEciIpJUUidRVJZz7O6XYN+meEciIpJUUidRHD7ASRt+Cx/8M96RiIgkldRJFPUns3XVk4hIe6ROoqi/PDZ1miwi0hlSZ6upbsZFRGKSOolCh55ERGKSOokiexhvTf0NnHRRvCMREUkqqZMo0jP5pPdx0HNAvCMREUkqqZMoDh8gd+ciKN0Y70hERJKKr4nCzGaa2QYz22xmtzZTfqWZvRd+LTWzib4FU76P0Vv+ALtX+bYKEZHuyLdEYWYB4G7gQmAscIWZjW1UbRtwjnNuAvBTYL5f8dSfzNblsSIi7eLnVnMasNk5t9U5VwksBGZHVnDOLXXOfRweXQbk+hZN3X0UIiLSLuk+Lns4sDNivBg4NUr9rwJ/a67AzOYB8wBycnIIBoPtDqZX+Q6mAWvWraN0f/vnT2ShUCimzySRqU3JQW1KDqFQqEPz+5komrthodmHQpjZDLxEcWZz5c65+YQPSxUWFrqioqL2R7N3HSyHcWPHwfgY5k9gwWCQmD6TBKY2JQe1KTl0NPH5eeipGBgRMZ4LlDSuZGYTgPuB2c65/b5FM2g0b0y/H074tG+rEBHpjvxMFMuBMWaWb2aZwOXAosgKZnYc8CRwlXPO3+tWAxkcyRoCmb19XY2ISHfj26En51y1md0I/B0IAA8459aY2fXh8nuBHwKDgN+Z1wdTtXOu0JeADu3h+O3/B/uGw+AxvqxCRKQ78vMcBc6554HnG027N2L4WuBaP2OoF/qQ/O1/htLPKFGIiLRD6txU4HQfhYhILFJnq1n/PAr1Hisi0h6pkyjUzbiISExSJ1HU5wklChGR9vD1ZHZCGTqRf5zxCGeOnBHvSKQVtbWO6lpHRVUNNeHhmsiXc9TUhN9ra6l13imoWueodQ4XHnc4asPTvWkOF15+bbjcNZ4Xr15trffbInLe2ohluvD0unI4+lvEhSfU310aHlhXXMWe5Tua1PfmqZvmGo03rOAa12+8rmbmravT1nU1Lo+c1tjWrZWscZubL2xm3a1pZ/Xm7+Dt4PK3b6/k3epN4eUnVvxT8wZw1pgh7Zyr41InUQTSqc7oA+mZ8Y4kodTUOj6prOZwVQ0VlbUcrqo5Ol5Vw+HKWj6prPaGq2qoqKqlsrqWqppajoTf68Yra2qprHZU1tRS1WCa915VU0tNjbfhr3XNJIBwEqj/sr34Qlw/G1+sfj/eEXS+TRviHUHn25yYjyO4/pxRShS+OrCDkVsWwPhcGDw63tF0KuccBw9XsT90hP3llewPHWFfqJKDh6soq6jiUEU1ZYe990MVVZTVvR/2EkJ7BdKMjICRGUgjMz2NjMj3QBoZ6Wn0CKTRIyON7Kx0MsLTMgNppKcZgcYvMwIB7z09zUhLM3Z8sJ3Ro0Z64+HpXv20+jp172kGaea9g/dudnQ64Xcjol7ktLSG86aZYXXz1s2TBla/7Lrle/N7q/DqNj6yWTduZix74w1OO+20htMjzpkdnUaDgcbLtvr61mg8Yr00rNzWeVtqR0sxv/baEs4++5ymlZuZr63ae3DY2rmC1moHlwQpOqcoYvn+xpMMUidRlJVw3M6n4MBVSZUoKqtr2VNWwa4Dhyk5cJjdB48Olx46wv5QJaWHKqj5+4vNzp8Z8DbWfXtmkJ2VTnZWOjl9s7xpWRn0yUqnd2Y6WZkBemYE6BV+z6obDo/3jJgeSPP/ixAMllBUlDx/p7YY1DONYf17xjuMTpWeZmSmd69Tnd4Pg+63se+I1EkUCX4fxb7QETZ8eIitpSG27itna2k5W/eF2PXxYWobHcgc2DuTYf2zyOmbxbhhfQnt38OUsaMZ1CeTQb171L/375VBj/S0bvkLR0S6TgolisS5j6L00BFWfvAxa0oOsqakjDUlB9lTdqS+vFdmgJFDejNpxAA+NzmX3AE9Gd6/J0P7ZTG0X096ZgYaLC8YDFJ01siuboaIpIjUSRRxvI9i76EK/rl5H29t+4g3t33E1tJyANIMRh/Th9NHDWbcsL6cPLQvo4b0IadvD+0FiEjCSJ1EUX/oyf8NsHOOtbvLeGXdXl5Zt4d3iw8CkJ2VzrS8gXyhcART8wYyblhfsjICrSxNRCS+UidR5J3JkrOf5Jzjm302UqfY+dEnPL1qF0++s4utpeWYwaQR/fnOp0/knBOGcPLQvl1yIlhEpDOlTqIww6UFIK1zT2Y753ht0z4e+Mc2lmwsBWBa/kCuO2skF4zNYXCfHp26PhGRrpY6iaJ0I2M23gunjIBBozq8uOqaWp58ZxfzX9vK5r0hjsnuwbfOP4FLC4aTO6BXJwQsIpIYUidRlBUzvORvUP7NDiWK2lrHs+/v5lcvbWTrvnLGD+/Lry6bxEWnDO1215OLiEAqJYq6y2M7cNXT+g/L+MFTq1nxwcecmJPN/KsKuGBsjq5QEpFuLYUSRfg9hhvujlTXcOdLG7n/9W30zUrn9ktP4fMFI3RiWkRSQgolithuuNu89xA3PbqKdbvL+EJhLt+78GQG9FbHgiKSOlInUVgaNWmZBNqxR/Hk28X8v6fep1dmOvdfXcj5Y3N8DFBEJDGlTqIYcz6vn/0YRcOntFrVOcevXt7EXa9s4rSRg7jr8kkc0zerC4IUEUk8qZMo2qiyupbvPfk+T7xdzOcLcvmvz52iq5lEJKWlTqIoWcVJ6/4XJubDgOObrXLwcBVf/9NKlm7Zz80XnMBN547WFU0ikvJSJ1EcLObYPUGoONhscfHHn/AvC5azbV85d35hIpdMye3a+EREElTqJIooVz2t3nWQaxYsp6Kqhoeumcbpowd3cXAiIonL14PvZjbTzDaY2WYzu7WZcjOzX4fL3zOz1s80x6z5bsYXvVvC3HvfIDOQxhNfP11JQkSkEd/2KMwsANwNXAAUA8vNbJFzbm1EtQuBMeHXqcA94ffO16ib8V0HDvPrlzfxfyt2Unj8AH73pSkck60rm0REGvPz0NM0YLNzbiuAmS0EZgORiWI28LBzzgHLzKy/mQ11zu3u9GgCmRy0bG58cCXrqnezL3SE9DTja+eM5N8vOFFXNomItMDPRDEc2BkxXkzTvYXm6gwHGiQKM5sHzAPIyckhGAzGEE4v/jTwXqpq0xkXqGHIsAxOPTadIT33sPQfe2JYXuIIhUIxfiaJS21KDmpTcgiFQh2a389E0dx1pS6GOjjn5gPzAQoLC11RUVGMIQWJfd7EFQx2v3apTclBbUoOHU18fh5vKQZGRIznAiUx1BERkTjyM1EsB8aYWb6ZZQKXA4sa1VkEXB2++mk6cNCX8xMiIhIz3w49OeeqzexG4O9AAHjAObfGzK4Pl98LPA9cBGwGPgGu8SseERGJja833DnnnsdLBpHT7o0YdsANfsYgIiIdo2tCRUQkKiUKERGJSolCRESiUqIQEZGozLkm97clNDMrBT6IcfbBwL5ODCdRdMd2qU3JQW1KDoOB3s65IbHMnHSJoiPMbIVzrjDecXS27tgutSk5qE3JoaNt0qEnERGJSolCRESiSrVEMT/eAfikO7ZLbUoOalNy6FCbUuochYiItF+q7VGIiEg7KVGIiEhUKZMozGymmW0ws81mdmu842krM3vAzPaa2eqIaQPN7CUz2xR+HxBR9r1wGzeY2afjE3V0ZjbCzBab2TozW2Nm/xaenrTtMrMsM3vLzN4Nt+kn4elJ26Y6ZhYws3fM7NnweFK3ycy2m9n7ZrbKzFaEpyV7m/qb2eNmtj78vTqtU9vknOv2L7xuzrcAI4FM4F1gbLzjamPsZwNTgNUR0+4Abg0P3wrcHh4eG25bDyA/3OZAvNvQTJuGAlPCw9nAxnDsSdsuvKc19gkPZwBvAtOTuU0RbbsZ+DPwbDf5/9sODG40Ldnb9BBwbXg4E+jfmW1KlT2KacBm59xW51wlsBCYHeeY2sQ59xrwUaPJs/H+MQi/z4mYvtA5d8Q5tw3vOR/TuiLO9nDO7XbOvR0ePgSsw3tWetK2y3nqHkycEX45krhNAGaWC3wGuD9iclK3qQVJ2yYz64v3g/IPAM65SufcATqxTamSKIYDOyPGi8PTklWOCz8JMPx+THh60rXTzPKAyXi/wJO6XeFDNKuAvcBLzrmkbxPwK+AWoDZiWrK3yQEvmtlKM5sXnpbMbRoJlAIPhg8R3m9mvenENqVKorBmpnXH64KTqp1m1gd4Avimc64sWtVmpiVcu5xzNc65SXjPfp9mZuOjVE/4NpnZZ4G9zrmVbZ2lmWkJ1aawM5xzU4ALgRvM7OwodZOhTel4h6fvcc5NBsrxDjW1pN1tSpVEUQyMiBjPBUriFEtn2GNmQwHC73vD05OmnWaWgZckHnHOPRmenPTtAgjv9geBmSR3m84ALjaz7XiHa881sz+R3G3COVcSft8LPIV32CWZ21QMFIf3YAEex0scndamVEkUy4ExZpZvZpnA5cCiOMfUEYuAL4eHvww8HTH9cjPrYWb5wBjgrTjEF5WZGd7x1HXOuTsjipK2XWY2xMz6h4d7AucD60niNjnnvuecy3XO5eF9Z151zn2JJG6TmfU2s+y6YeBTwGqSuE3OuQ+BnWZ2YnjSecBaOrNN8T5b34VXBVyEd3XNFuD78Y6nHXE/CuwGqvB+CXwVGAS8AmwKvw+MqP/9cBs3ABfGO/4W2nQm3q7ue8Cq8OuiZG4XMAF4J9ym1cAPw9OTtk2N2lfE0auekrZNeMfz3w2/1tRtC5K5TeEYJwErwv9/fwUGdGab1IWHiIhElSqHnkREJEZKFCIiEpUShYiIRKVEISIiUSlRiIhIVEoUktLCvW5+I2J8mJk97tO65pjZD1soC4Xfh5jZC36sXyRWShSS6voD9YnCOVfinPu8T+u6BfhdtArOuVJgt5md4VMMIu2mRCGp7ufAqPCzCX5hZnkWfvaHmX3FzP5qZs+Y2TYzu9HMbg53vLbMzAaG640ysxfCncy9bmYnNV6JmZ0AHHHO7QuP55vZG2a23Mx+2qj6X4ErfW21SDsoUUiquxXY4pyb5Jz7TjPl44Ev4vUH9DPgE+d1vPYGcHW4znzgJudcAfBtmt9rOAN4O2L8LrxO3KYCHzaquwI4K8b2iHS69HgHIJLgFjvvmRmHzOwg8Ex4+vvAhHAPuKcDj3ldWAHeA2EaG4rXFXSdM4BLw8N/BG6PKNsLDOuc8EU6TolCJLojEcO1EeO1eN+fNOCA87oXj+Yw0K/RtJb6z8kK1xdJCDr0JKnuEN7jWGPivOdobDOzueD1jGtmE5upug4YHTH+T7weWaHp+YgT8DoWFEkIShSS0pxz+4F/mtlqM/tFjIu5EviqmdX1SNrcY3ZfAybb0eNT/4b30JzlNN3TmAE8F2MsIp1OvceKdBEzuwt4xjn3civ1XgNmO+c+7prIRKLTHoVI1/kvoFe0CmY2BLhTSUISifYoREQkKu1RiIhIVEoUIiISlRKFiIhEpUQhIiJRKVGIiEhU/x83qy1trjG0nAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.plot(t[1:], step, label='analytic')\n", "plt.plot(t, stepnum, '--', label='numerical')\n", "plt.xlabel('time (d)')\n", "plt.ylabel('step (m)')\n", "plt.grid()\n", "plt.legend();" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.9.12" } }, "nbformat": 4, "nbformat_minor": 5 }