lectures/waveEquation1D.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "59dcf65a",
   "metadata": {},
   "source": [
    "# The wave equation\n",
    "The wave equation describes motion of a string.\n",
    "![string.png](string.png)\n",
    "The wave equation derives from Newton's law and Hook's law\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "\\rho u_{tt}(x,t) & = & \\sigma_x(x,t) + f(x,t), \\label{eq:1} \\\\\n",
    "\\sigma(x,t)           & = & \\kappa u_x (x,t)\n",
    "            \\label{eq:2}\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "Here $u(x,t)$ is the displacement (from equilibrium) of the string\n",
    "and $\\sigma$ is the stress (force per length unit) in the string.\n",
    "$\\kappa(x)$ is the stiffness and $\\rho$ is the density.\n",
    "$f(x,t)$ is a driving force (source).\n",
    "\n",
    "Differentiating the second equation above  with respect to $x$ and inserting into\n",
    "the first equation we get\n",
    "\\begin{eqnarray}\n",
    "  \\frac{\\rho}{\\kappa} u_{tt} = u_{xx} + \\frac{f}{\\kappa}.\n",
    "\\end{eqnarray}\n",
    "Where we assume that $\\rho$ is constant. The\n",
    "wave velocity is given by $c^{-2} = \\rho/\\kappa$, so we have\n",
    "\\begin{eqnarray}\n",
    "  c^{-2}(x) u_{tt}(x,t) = u_{xx}(x,t) + s(x,t).\n",
    "\\end{eqnarray}\n",
    "where $s(x,t) = \\frac{f(x,t)}{\\kappa}$.\n",
    "\n",
    "## Numerical Grid\n",
    "We create a grid out of the $x$ and $t$ coordinates as follows\n",
    "\\begin{eqnarray}\n",
    "  x_i & = & i\\Delta x = i h, \\, i=1,\\cdots,n\\\\\n",
    "  t_j & = & j\\Delta t = j k, \\, j=1,\\cdots,m.\n",
    "\\end{eqnarray}\n",
    "We want to compute an approximate solution, $u_{i,j}$ to the true solution\n",
    "$u(x_i,t_j)$ for all gridpoints in the grid.\n",
    "\n",
    "![grid.png](grid.png)\n",
    "\n",
    "## Finite Difference Solution\n",
    "The second order derivatives are approximated with a central difference\n",
    "\\begin{eqnarray}\n",
    "\\frac{u_{i,j+1} -2 u_{i,j} + u_{i,j-1}}{k^2} =\n",
    "c^2_i\\frac{u_{i+1,j} - 2u_{i,j} + u_{i-1,j}}{h^2} + s_{i,j}\n",
    "\\end{eqnarray}\n",
    "This is equal to\n",
    "\\begin{eqnarray}\n",
    "u_{i,j+1} -2 u_{i,j} + u_{i,j-1} =\n",
    "\\frac{k^2c^2_i}{h^2}(u_{i+1,j} - 2u_{i,j} + u_{i-1,j}) + s_{i,j}k^2\n",
    "\\end{eqnarray}\n",
    "Introduce $r^2_i = \\left(\\frac{kc_i}{h}\\right)^2$ to get\n",
    "\\begin{eqnarray}\n",
    "u_{i,j+1} -2 u_{i,j} + u_{i,j-1} =\n",
    "r^2_i(u_{i+1},j - 2u_{i,j} + u_{i-1,j}) + s_{i,j} k^2\n",
    "                                        \\label{eq:diff}\n",
    "\\end{eqnarray}\n",
    "\n",
    "which is solved by\n",
    "\\begin{eqnarray}\n",
    "u_{i,j+1} = 2 u_{i,j} - u_{i,j-1}\n",
    "+r^2_i(u_{i+1},j - 2u_{i,j} + u_{i-1,j}) + s_{i,j} k^2\n",
    "                                        \\label{eq:sol}\n",
    "\\end{eqnarray}\n",
    "\n",
    "## Stability\n",
    "Stability condition\n",
    "\\begin{eqnarray}\n",
    " r = ck/h \\le 1\n",
    "\\end{eqnarray}\n",
    "\n",
    "If $r=1$ the numerical solution is exact!\n",
    "An exact solution to the wave equation is a \n",
    "plane monochromatic wave with frequency $\\omega$ and\n",
    "wavenumber $k_x$ which are constants.\n",
    "\\begin{eqnarray}\n",
    "u(x,t) = \\exp i(k_x x - \\omega t)\n",
    "\\end{eqnarray}\n",
    "Prove that this is a solution!\n",
    "\n",
    "We will now find the stability condition by using\n",
    "the plane wave as a trial solution in the\n",
    "finite-difference solution.\n",
    "\n",
    "First, to avoid confusion, we replace $i$ and $j$ in the\n",
    "finite-difference solution above with indices $m$ (x) and \n",
    "$n$ (time). We then have the displacement $u$ as:\n",
    "    \n",
    "\\begin{eqnarray}\n",
    "u_{m,n} = \\exp i(k_x m h - \\omega n k)\n",
    "\\end{eqnarray}\n",
    "\n",
    "If we insert the above expression for $u$ into the \n",
    "finite-difference solution above (changing $i$ to $m$ \n",
    "and $j$ to $n$ one gets:\n",
    "                                  \n",
    "\\begin{eqnarray}\n",
    "  \\exp(-i\\omega k )  = 2 -\\exp(i\\omega k)   \n",
    "                    + r^2[\\exp(ik_x h)- 2 + \\exp(-ik_x h)]\n",
    "\\end{eqnarray}\n",
    "                                  \n",
    "This can be rearranged to:\n",
    "\\begin{eqnarray}\n",
    "\\exp(-i\\omega k) + \\exp(i\\omega k) = (2-2r^2) +\n",
    "                                     r^2[\\exp(ik_x m h)] + \\exp(-ik_x m h)]\n",
    "\\end{eqnarray}\n",
    "                                  \n",
    "Since $2\\cos(x) = \\exp(ix)+\\exp(-ix)$\n",
    "\n",
    "the relation above can be rewritten to\n",
    "                                  \n",
    "\\begin{eqnarray}\n",
    "\\cos(\\omega k) = (1-r^2) + r^2\\cos(k_x m h) \n",
    "\\end{eqnarray}\n",
    "                \n",
    "                           \n",
    "1. For $r=1$ this gives $\\omega h = k_x m $ which\n",
    "   is the so-called dispersion relation for a plane \n",
    "   wave, which means that $r=1$ corresponds to the exact solution.\n",
    "                            \n",
    "2. For $r < 1$ the finite-difference solution is stable, but not exact.\n",
    "\n",
    "3. For $r>1$ the finite-difference solution is unstable.\n",
    "                                  \n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "1451d2bd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAGwCAYAAABFFQqPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/av/WaAAAACXBIWXMAAA9hAAAPYQGoP6dpAABA/UlEQVR4nO3de3RU9b3//9dMLpNESAKGJIDRgCioCAh8SaNSrURBPShd1lJKERGxcsALoa2kVfByNNgKUimaI4qXs1RovR1+gmk5gUipUSRIEQVURKFIAjTmCuQys39/JLMz12QGctvm+VhrFpm998x8Zsc2r/X+3GyGYRgCAADohuyd3QAAAIDOQhACAADdFkEIAAB0WwQhAADQbRGEAABAt0UQAgAA3RZBCAAAdFuRnd2Ars7lcunbb79Vz549ZbPZOrs5AAAgBIZhqKqqSv369ZPdHrzuQxBqxbfffqu0tLTObgYAADgFBw8e1FlnnRX0PEGoFT179pTUeCPj4+M7uTUAACAUlZWVSktLM/+OB0MQaoW7Oyw+Pp4gBACAxbQ2rIXB0gAAoNsiCAEAgG6LIAQAALotghAAAOi2CEIAAKDbIggBAIBuiyAEAAC6LYIQAADotghCAACg2yIIAQCAbstSQWjz5s2aOHGi+vXrJ5vNprfffrvV1xQWFmrkyJFyOBwaNGiQXnzxxXZvJwAAsAZLBaGamhoNHz5cK1asCOn6/fv36/rrr9ePfvQj7dixQ/fee69uv/12/fWvf23nlgIAACuw1Kar1157ra699tqQr8/Ly9OAAQO0ZMkSSdIFF1ygLVu26Mknn9T48ePbq5kdzukyFGFveVM5AADgz1IVoXAVFRUpKyvL69j48eNVVFQU9DW1tbWqrKz0enRlee/t0+D739X2A991dlMAALCc73UQKikpUUpKitexlJQUVVZW6sSJEwFfk5ubq4SEBPORlpbWEU09ZYvf3aMGl6EFb+zs7KYAAGA53+sgdCpycnJUUVFhPg4ePNjZTQpJdCS/SgAAwmWpMULhSk1NVWlpqdex0tJSxcfHKzY2NuBrHA6HHA5HRzSvTTU4jc5uAgAAlvO9LiNkZmaqoKDA69iGDRuUmZnZSS1qP/VOV2c3AQAAy7FUEKqurtaOHTu0Y8cOSY3T43fs2KEDBw5IauzWuuWWW8zr77zzTn311Vf6zW9+oz179ujpp5/Wn//8Z82bN68zmg8AALoYSwWhbdu26ZJLLtEll1wiScrOztYll1yihQsXSpIOHz5shiJJGjBggNatW6cNGzZo+PDhWrJkiZ577rnv1dR5N5uN6fMAAITLUmOErrzyShlG8LEwgVaNvvLKK/Xxxx+3Y6s6j9PVfC+IQQAAhM9SFSF4q2toHhdEQQgAgPARhCystsFp/myjJgQAQNgIQhbmWRFyttBlCAAAAiMIWVitRxDyDEUAACA0BCEL8+waIwgBABA+gpCFeVaEWFARAIDwEYQsjK4xAABOD0HIwjzXEWpwMVgaAIBwEYQszHOjVSdBCACAsBGELMy7IkTXGAAA4SIIWZjn2kEuQ3JRFQIAICwEIQtz+lSBWFQRAIDwEIQszHOMkMQ4IQAAwkUQsjDf4MPMMQAAwkMQsjDf4ON0EoQAAAgHQcjCXIZvRYiZYwAAhIMgZGGMEQIA4PQQhCyMMUIAAJwegpCF+Y0RIggBABAWgpCF+a4jREUIAIDwEIQszLcC5BuMAABAywhCFuZbAaIiBABAeAhCFuY3WJp1hAAACAtByMIYLA0AwOkhCFmY727zdI0BABAegpCFURECAOD0EIQszH9BRWaNAQAQDoKQhVERAgDg9BCELIwFFQEAOD0EIQtz+vSEMX0eAIDwEIQszLciRNcYAADhIQhZmG9XmGEQhAAACAdByML89hojCAEAEBaCkIUxawwAgNNDELIw35WlKQgBABAegpCFURECAOD0WC4IrVixQunp6YqJiVFGRoa2bt3a4vXLli3T4MGDFRsbq7S0NM2bN08nT57soNa2L8YIAQBweiwVhNasWaPs7GwtWrRI27dv1/DhwzV+/HgdOXIk4PWvvvqqFixYoEWLFmn37t16/vnntWbNGv32t7/t4Ja3D98gxKwxAADCY6kgtHTpUs2aNUszZszQhRdeqLy8PMXFxWnVqlUBr3///fd12WWX6ec//7nS09N1zTXXaMqUKS1WkWpra1VZWen16Kp8K0C+CywCAICWWSYI1dXVqbi4WFlZWeYxu92urKwsFRUVBXzNpZdequLiYjP4fPXVV1q/fr2uu+66oJ+Tm5urhIQE85GWlta2X6QN+VaA6BoDACA8kZ3dgFAdO3ZMTqdTKSkpXsdTUlK0Z8+egK/5+c9/rmPHjunyyy+XYRhqaGjQnXfe2WLXWE5OjrKzs83nlZWVXTYM+XaN+c4iAwAALbNMRehUFBYW6rHHHtPTTz+t7du3680339S6dev0yCOPBH2Nw+FQfHy816Or8s09LipCAACExTIVoaSkJEVERKi0tNTreGlpqVJTUwO+5oEHHtC0adN0++23S5Iuvvhi1dTU6I477tDvfvc72e3WzoG+wYfp8wAAhMcySSA6OlqjRo1SQUGBeczlcqmgoECZmZkBX3P8+HG/sBMRESHp+zHDyh2E7Dbv5wAAIDSWqQhJUnZ2tqZPn67Ro0drzJgxWrZsmWpqajRjxgxJ0i233KL+/fsrNzdXkjRx4kQtXbpUl1xyiTIyMvTll1/qgQce0MSJE81AZGXuClBUhF21DS6/rjIAANAySwWhyZMn6+jRo1q4cKFKSko0YsQI5efnmwOoDxw44FUBuv/++2Wz2XT//ffr0KFD6tOnjyZOnKhHH320s75Cm3IHH3cQomsMAIDw2IzvQx9RO6qsrFRCQoIqKiq63MDpnzzzvrZ9850S46JUfrxe868+X3eNO6+zmwUAQKcL9e+3ZcYIwZ97TFBkUxWMghAAAOEhCFlYc9dY42hpFlQEACA8BCELMytCTUGIBRUBAAgPQcjC3EEoqqlrjIoQAADhIQhZmHuT1agI9xghghAAAOEgCFmYQdcYAACnhSBkYe51gyKbKkLuChEAAAgNQcjCmscI2byeAwCA0BCELMxzZenG5wQhAADCQRCyMN/p82yxAQBAeAhCFmZ2jUWwsjQAAKeCIGRhrqbB0ZF2Zo0BAHAqCEIW5lsRYkFFAADCQxCysObp81SEAAA4FQQhC3Pnnubd5wlCAACEgyBkYe7gEx3p3n2+M1sDAID1EIQszJw+764I0TUGAEBYCEIW5jdGiK4xAADCQhCyMMNnZWkWVAQAIDwEIQtr7hqjIgQAwKkgCFmY7+7zFIQAAAgPQcjC3AWgaPYaAwDglBCELMxp+FaECEIAAISDIGRhvmOEqAgBABAegpBFGYbhN2uMihAAAOEhCFmUZ/Gnea+xTmoMAAAWRRCyKM9usCg7u88DAHAqCEIW5dkNxsrSAACcGoKQRXlmHnOMEIOlAQAIC0HIojy7waLc6whREQIAICwEIYvy6hozd5/vrNYAAGBNBCGL8uwGY4wQAACnhiBkUZ7DgaLZfR4AgFNCELIoz9ATYWeMEAAAp4IgZFFGU+ix25qDEDkIAIDwWC4IrVixQunp6YqJiVFGRoa2bt3a4vXl5eWaM2eO+vbtK4fDofPPP1/r16/voNa2H3f1J8Juk83GGCEAAE5FZGc3IBxr1qxRdna28vLylJGRoWXLlmn8+PHau3evkpOT/a6vq6vT1VdfreTkZL3++uvq37+/vvnmGyUmJnZ849uYu2fMZrOpqSDEGCEAAMJkqSC0dOlSzZo1SzNmzJAk5eXlad26dVq1apUWLFjgd/2qVatUVlam999/X1FRUZKk9PT0jmxyu3HPGouw2egaAwDgFFmma6yurk7FxcXKysoyj9ntdmVlZamoqCjga9auXavMzEzNmTNHKSkpGjp0qB577DE5nc6gn1NbW6vKykqvR1fk8hgjZG/qGqMiBABAeCwThI4dOyan06mUlBSv4ykpKSopKQn4mq+++kqvv/66nE6n1q9frwceeEBLlizRf/3XfwX9nNzcXCUkJJiPtLS0Nv0ebcWdeew2mxmEGCMEAEB4LBOEToXL5VJycrKeffZZjRo1SpMnT9bvfvc75eXlBX1NTk6OKioqzMfBgwc7sMWhc1d/7HabmhaWJggBABAmy4wRSkpKUkREhEpLS72Ol5aWKjU1NeBr+vbtq6ioKEVERJjHLrjgApWUlKiurk7R0dF+r3E4HHI4HG3b+HbgNX3erAh1ZosAALAey1SEoqOjNWrUKBUUFJjHXC6XCgoKlJmZGfA1l112mb788ku5PDbh+vzzz9W3b9+AIchKAk2fZ4wQAADhsUwQkqTs7GytXLlSL730knbv3q3Zs2erpqbGnEV2yy23KCcnx7x+9uzZKisr0z333KPPP/9c69at02OPPaY5c+Z01ldoM+5sZ/OYNUbXGAAA4bFM15gkTZ48WUePHtXChQtVUlKiESNGKD8/3xxAfeDAAdntzdkuLS1Nf/3rXzVv3jwNGzZM/fv31z333KP77ruvs75Cm3GHngiPdYRcVIQAAAiLpYKQJM2dO1dz584NeK6wsNDvWGZmpj744IN2blXHCzR9nhwEAEB4LNU1hmbm9Hm7TXa6xgAAOCUEIYsyp897do0RhAAACAtByKIMj1ljTJ8HAODUEIQsyl0RstnE9HkAAE4RQciiPLfYcE+fl5orRQAAoHUEIYsKNH1eoioEAEA4CEIW5Q5CNpvMWWONxzurRQAAWA9ByKLclZ8Ie/Pu8xIzxwAACAdByKIMzzFCBCEAAE4JQciizJWl7TbZGCMEAMApIQhZVPOCivLpGuusFgEAYD0EIYtyB54Ips8DAHDKCEIW1bzpKtPnAQA4VQQhi2oeI9S4srTN3G+sExsFAIDFEIQsynPTVc9/mTUGAEDoCEIW5Tl9XpLHxqsEIQAAQkUQsiizItQ0QMjdNcYYIQAAQkcQsqjmwdKNz90zxygIAQAQOoKQRRke0+el5i4yKkIAAISOIGRRTnPTVXcQajzOGCEAAEJHELIod+CJaPoNuscKURACACB0BCGLcjF9HgCA00YQsih35cddCSIIAQAQPoKQRfkvqOh9HAAAtI4gZFHmGCGmzwMAcMoIQhbluemq579UhAAACB1ByKL8xgjZ3ccJQgAAhIogZFHNY4TU9C+DpQEACBdByKIMn66x5k1XO61JAABYDkHIony7xth0FQCA8BGELCrYpqt0jQEAEDqCkEUFW1maHAQAQOgIQhZldo3Z3F1jTJ8HACBcBCGL8l1HKILp8wAAhI0gZFHNFSH3v4wRAgAgXAQhizIrQr6brro6rUkAAFiO5YLQihUrlJ6erpiYGGVkZGjr1q0hvW716tWy2WyaNGlS+zawg/gPlm487qQiBABAyE4pCJWXl+u5555TTk6OysrKJEnbt2/XoUOH2rRxvtasWaPs7GwtWrRI27dv1/DhwzV+/HgdOXKkxdd9/fXX+tWvfqWxY8e2a/s6km/XWPOmqwQhAABCFXYQ2rlzp84//3w9/vjjeuKJJ1ReXi5JevPNN5WTk9PW7fOydOlSzZo1SzNmzNCFF16ovLw8xcXFadWqVUFf43Q6NXXqVD300EMaOHBgq59RW1uryspKr0dX5DtYunnWWKc1CQAAywk7CGVnZ+vWW2/VF198oZiYGPP4ddddp82bN7dp4zzV1dWpuLhYWVlZ5jG73a6srCwVFRUFfd3DDz+s5ORkzZw5M6TPyc3NVUJCgvlIS0s77ba3B78FFRksDQBA2MIOQh999JF++ctf+h3v37+/SkpK2qRRgRw7dkxOp1MpKSlex1NSUoJ+7pYtW/T8889r5cqVIX9OTk6OKioqzMfBgwdPq93txW+wNNPnAQAIW2S4L3A4HAG7iz7//HP16dOnTRrVFqqqqjRt2jStXLlSSUlJIb/O4XDI4XC0Y8vahu+CikyfBwAgfGEHoRtuuEEPP/yw/vznP0tqHJty4MAB3XfffbrpppvavIFuSUlJioiIUGlpqdfx0tJSpaam+l2/b98+ff3115o4caJ5zNU0tzwyMlJ79+7Vueee227tbW/Ns8bU9C/T5wEACFfYXWNLlixRdXW1kpOTdeLECV1xxRUaNGiQevbsqUcffbQ92ihJio6O1qhRo1RQUGAec7lcKigoUGZmpt/1Q4YM0SeffKIdO3aYjxtuuEE/+tGPtGPHji479idU7sqPjenzAACcsrArQgkJCdqwYYO2bNminTt3qrq6WiNHjvQaxNxesrOzNX36dI0ePVpjxozRsmXLVFNToxkzZkiSbrnlFvXv31+5ubmKiYnR0KFDvV6fmJgoSX7HrcjdNeaeNs/0eQAAwhd2EHK7/PLLdfnll7dlW1o1efJkHT16VAsXLlRJSYlGjBih/Px8cwD1gQMHZLdbbo3IU+I7a4zp8wAAhC+kIPTUU0+F/IZ33333KTcmFHPnztXcuXMDnissLGzxtS+++GLbN6iT+K4szfR5AADCF1IQevLJJ72eHz16VMePHze7msrLyxUXF6fk5OR2D0Jo5O4aM8cIMX0eAICwhdSPtH//fvPx6KOPasSIEdq9e7fKyspUVlam3bt3a+TIkXrkkUfau71o4g48EX6zxghCAACEKuwBNQ888ICWL1+uwYMHm8cGDx6sJ598Uvfff3+bNg7BuQs/vrvPO8lBAACELOwgdPjwYTU0NPgddzqdfmv8oP04Xd7T55k1BgBA+MIOQuPGjdMvf/lLbd++3TxWXFys2bNnd8gUejRq7hpzb7raeNxJ1xgAACELOwitWrVKqampGj16tLkdxZgxY5SSkqLnnnuuPdqIAJq32HD/a/M6DgAAWhf2OkJ9+vTR+vXr9fnnn2vPnj2SGldxPv/889u8cQjOMJg+DwDA6TrlBRXPP/98wk8ncppbbDQ+N6fPUxICACBkYQeh2267rcXzq1atOuXGIHS+W2zQNQYAQPjCDkLfffed1/P6+nrt2rVL5eXluuqqq9qsYWiZb9dY8/R5khAAAKEKOwi99dZbfsdcLpdmz56tc889t00ahda5fLrGmD4PAED42mSHUrvdruzsbL+tONB+nD57jTF9HgCA8LXZVu379u0LuNAi2ofvGKEIxggBABC2sLvGsrOzvZ4bhqHDhw9r3bp1mj59eps1DC1rHiPU+Ny91QbT5wEACF3YQejjjz/2em6329WnTx8tWbKk1RllaDt+u8+z6SoAAGELOwht2rSpPdqBMPmOEXJXhshBAACELuwxQldddZXKy8v9jldWVjJ9vgO5u8Yimn6DdlaWBgAgbGEHocLCQtXV1fkdP3nypP7+97+3SaPQOr+uMcYIAQAQtpC7xnbu3Gn+/Nlnn6mkpMR87nQ6lZ+fr/79+7dt6xCUy29BxcbjTJ8HACB0IQehESNGyGazyWazBewCi42N1fLly9u0cQiueYxQ43OmzwMAEL6Qg9D+/ftlGIYGDhyorVu3qk+fPua56OhoJScnKyIiol0aCX/uHrAI364xkhAAACELOQidc845khq300Dna95iw3fTVYIQAAChCikIrV27Vtdee62ioqK0du3aFq+94YYb2qRhaJnLd0FF9xghghAAACELKQhNmjRJJSUlSk5O1qRJk4JeZ7PZ5HQ626ptaIG7B8zdJda86WpntQgAAOsJKQh5dofRNdY1+M4ac3eRMWsMAIDQtdmmq+hYvl1jETbv4wAAoHUhVYSeeuqpkN/w7rvvPuXGIHTuwpzdZ9YYOQgAgNCFFISefPLJkN7MZrMRhDoIXWMAAJy+kILQ/v3727sdCJMZhJo6NyOYPg8AQNhOa4yQYRjm5p/oWOasMb/d5/l9AAAQqlMKQs8//7yGDh2qmJgYxcTEaOjQoXruuefaum1ogcvls9eYnS02AAAIV8grS7stXLhQS5cu1V133aXMzExJUlFRkebNm6cDBw7o4YcfbvNGwp//goqMEQIAIFxhB6FnnnlGK1eu1JQpU8xjN9xwg4YNG6a77rqLINRB/BdUdB8nCAEAEKqwu8bq6+s1evRov+OjRo1SQ0NDmzQKrfOdNcZeYwAAhC/sIDRt2jQ988wzfsefffZZTZ06tU0ahdY1jxFS07/u3ec7q0UAAFjPaQ2Wvv3223X77bfr4osv1sqVK2W325WdnW0+2sOKFSuUnp6umJgYZWRkaOvWrUGvXblypcaOHatevXqpV69eysrKavF6K/GfNdY0RoiKEAAAIQt7jNCuXbs0cuRISdK+ffskSUlJSUpKStKuXbvM69wL/LWlNWvWKDs7W3l5ecrIyNCyZcs0fvx47d27V8nJyX7XFxYWasqUKbr00ksVExOjxx9/XNdcc40+/fRT9e/fv83b15Ga1xHynj7PcgYAAIQu7CC0adOm9mhHSJYuXapZs2ZpxowZkqS8vDytW7dOq1at0oIFC/yuf+WVV7yeP/fcc3rjjTdUUFCgW265pUPa3F4MsyLU9C/T5wEACJtlNl2tq6tTcXGxsrKyzGN2u11ZWVkqKioK6T2OHz+u+vp69e7dO+g1tbW1qqys9Hp0Rc4gg6WZPg8AQOjCrgidPHlSy5cv16ZNm3TkyBG5fEbnbt++vc0a5+nYsWNyOp1KSUnxOp6SkqI9e/aE9B733Xef+vXr5xWmfOXm5uqhhx46rbZ2BN9ZY+7p83SNAQAQurCD0MyZM/W3v/1NP/nJTzRmzJh2GQvUHhYvXqzVq1ersLBQMTExQa/LycnxGuhdWVmptLS0jmhiyBq3Nmn82d01ZmOwNAAAYQs7CL3zzjtav369LrvssvZoT1BJSUmKiIhQaWmp1/HS0lKlpqa2+NonnnhCixcv1v/93/9p2LBhLV7rcDjkcDhOu73tyTPrmBUhps8DABC2sMcI9e/fXz179myPtrQoOjpao0aNUkFBgXnM5XKpoKDA3OojkN///vd65JFHlJ+fH3AhSCvyrPqwoCIAAKcu7CC0ZMkS3Xffffrmm2/aoz0tys7O1sqVK/XSSy9p9+7dmj17tmpqasxZZLfccotycnLM6x9//HE98MADWrVqldLT01VSUqKSkhJVV1d3eNvbkmfYsdu9/yUIAQAQurC7xkaPHq2TJ09q4MCBiouLU1RUlNf5srKyNmucr8mTJ+vo0aNauHChSkpKNGLECOXn55sDqA8cOCC7vTnbPfPMM6qrq9NPfvITr/dZtGiRHnzwwXZrZ3sL1DXGrDEAAMIXdhCaMmWKDh06pMcee0wpKSkdPlh67ty5mjt3bsBzhYWFXs+//vrr9m9QJ3AF6BqLaBo1TUEIAIDQhR2E3n//fRUVFWn48OHt0R6EwLPqYzP3Gms6RxICACBkYY8RGjJkiE6cONEebUGIPHu/3JUgG4OlAQAIW9hBaPHixZo/f74KCwv173//2xKrMH/fGIG6xpg+DwBA2MLuGpswYYIkady4cV7HDcOQzWaT0+lsm5YhKJfXYGn3v1SEAAAIl6U2XUUj7zFCTbPGmD4PAEDYwg5CV1xxRdBzu3btOq3GIDTurjH3+CDJc/p8pzQJAABLOu3d56uqqvTss89qzJgxzCTrIC6ffcYkz+nzVIQAAAjVKQehzZs3a/r06erbt6+eeOIJXXXVVfrggw/asm0Iwt395bmGE9PnAQAIX1hdYyUlJXrxxRf1/PPPq7KyUj/96U9VW1urt99+WxdeeGF7tRE+3GOEPCtC5mBpVpYGACBkIVeEJk6cqMGDB2vnzp1atmyZvv32Wy1fvrw924Yg3EWfCJv/GCFyEAAAoQu5IvTuu+/q7rvv1uzZs3Xeeee1Z5vQCnfXmN0jCLnHCDFrDACA0IVcEdqyZYuqqqo0atQoZWRk6E9/+pOOHTvWnm1DEM1jhJqPuX9m01UAAEIXchD6wQ9+oJUrV+rw4cP65S9/qdWrV6tfv35yuVzasGGDqqqq2rOd8GBWhAJMn6cgBABA6MKeNXbGGWfotttu05YtW/TJJ59o/vz5Wrx4sZKTk3XDDTe0RxvhwxVgjBBdYwAAhO+01hEaPHiwfv/73+tf//qXXnvttbZqE1oRaPq8jenzAACE7bQXVJSkiIgITZo0SWvXrm2Lt0Mr3Burei2o6NE1xqKKAACEpk2CEDpWoFljnj8zXhoAgNAQhCzIFWivMbvN7zwAAGgZQciC3BUfm9fK0s0/M4UeAIDQEIQsqKUFFSWm0AMAECqCkAW5WthrTGLmGAAAoSIIWZC75yvQgoqN5wlCAACEgiBkQYFnjTWfN1wd3SIAAKyJIGRBdI0BANA2CEIWZHaN2Zg+DwDA6SAIWVCgrrHG503nmT4PAEBICEIW1Lz7vPfx5o1XO7pFAABYE0HIgoJVhNybsDJGCACA0BCELMi96arNJwi5N16lawwAgNAQhCzI3GvMOwc1jxGiIgQAQEgIQhYUaNaY1DxzjIIQAAChIQhZUPBZY01jhEhCAACEhCBkQcFmjbm7xgy6xgAACAlByIKCdY0xfR4AgPAQhCzIaG36PEkIAICQEIQsyB10fHJQ8/R5usYAAAiJ5YLQihUrlJ6erpiYGGVkZGjr1q0tXv+Xv/xFQ4YMUUxMjC6++GKtX7++g1raftwFnwh7kC02CEIAAITEUkFozZo1ys7O1qJFi7R9+3YNHz5c48eP15EjRwJe//7772vKlCmaOXOmPv74Y02aNEmTJk3Srl27OrjlbSvorDHGCAEAEJbIzm5AOJYuXapZs2ZpxowZkqS8vDytW7dOq1at0oIFC/yu/+Mf/6gJEybo17/+tSTpkUce0YYNG/SnP/1JeXl5AT+jtrZWtbW15vPKysp2+CbSov/dpe0HyhUZYVOU3a7ICJsi7DZFR9g1KKWHfpFxjtJ6xwV8bfMYIe/joUyf37q/TG8U/0vflNWorsGlBpchCkgAgM700I0XaeTZvTrlsy0ThOrq6lRcXKycnBzzmN1uV1ZWloqKigK+pqioSNnZ2V7Hxo8fr7fffjvo5+Tm5uqhhx5qkza35KtjNfrkUEXAcwV7jujVDw/o9Tsv1eDUnn7nncG22GhKRsGmz/9520H95vWdp9FqAADaXvXJhk77bMsEoWPHjsnpdColJcXreEpKivbs2RPwNSUlJQGvLykpCfo5OTk5XuGpsrJSaWlpp9HywH49frBuu2yA6p2NVZl6p0sNTkMnG5x6besB7TpUqUVrd2n1HZl+r23eYsN31ljjv4EqQt/V1Onh/+8zSdINw/tp3AXJio2KUGSETTbZ/K4HAKCjXNQvvtM+2zJBqKM4HA45HI52/5xhZyUGPXfVkGSNfXyTPviqTPuOVuvcPj28zhtBFlRsnjXm/57v7PxW1bUNGpLaU8smjzDHEwEA0J1ZZrB0UlKSIiIiVFpa6nW8tLRUqampAV+Tmpoa1vVdRd+EWP1g4JmSpMK9R/3Ou4OOb9eYvYXp8wV7GgeU3ziiPyEIAIAmlglC0dHRGjVqlAoKCsxjLpdLBQUFysz07z6SpMzMTK/rJWnDhg1Br+9KLhuUJEnafuA7v3Puri//BRUb//UNQoZhqPibxvcZe15SWzcVAADLslTXWHZ2tqZPn67Ro0drzJgxWrZsmWpqasxZZLfccov69++v3NxcSdI999yjK664QkuWLNH111+v1atXa9u2bXr22Wc782uExN1fuvuw/6y15jFC3seDbbFxuOKkqk42KMJu03kp3t1sAAB0Z5YKQpMnT9bRo0e1cOFClZSUaMSIEcrPzzcHRB84cEB2j4Ezl156qV599VXdf//9+u1vf6vzzjtPb7/9toYOHdpZXyFkF/RtDEL7j9XoZL1TMVER5jkjyF5jZteYTxLaW1IlSTq3zxlyREYIAAA0slQQkqS5c+dq7ty5Ac8VFhb6Hbv55pt18803t3Or2l5Sj2jFRUfoeJ1T35af0ECPAdPuipDfGCF74DFC+4/VSJIGJVMNAgDAk2XGCHU3NptNZ/WKlSQdKj/hdc4ZdEHFpvM+FaGSypOSpH4Jse3QUgAArIsg1IWd1atxZel/fecdhIJ1jQWbPv9tU5BKTYhph1YCAGBdBKEurG9TcDlccdLruDlrzG/T1cBdYyVNr+9LRQgAAC8EoS7szB6NCzuW1dR6HXcHoQif3557nLhvEHIHqb6JVIQAAPBEEOrCzjwjWpL07+o6r+PBttgItunqserGINWnR/uvmA0AgJUQhLqw3u4gVOMdhIJ1jTVvutp87GS9U7UNjbu09mp6PwAA0Igg1IU1V4R8usaCbrrqXxEqP14vSYq023RGNGsIAQDgiSDUhTWPEfLpGjPHCPl2jTWd9ygJfXe88bWJcVF+6w4BANDdEYS6MHfX2HfH69XgdJnH3T/6dY3Z/LvG3BWhhNiodmwpAADWRBDqwuJjmxf+rql1mj8HGyxtdo15JKGKE+6KEOODAADwRRDqwhyREYqObPwVVZ6sN48HHyzd+K9n15i7IpRIRQgAAD8EoS4uPqaxKlR1ssE8FmywdKBNV8tPNHWNxRGEAADwRRDq4nrGNAaYKo+KkCvogor+W2yYg6Vj6RoDAMAXQaiL69lUEaqu9agItbLFhuf0+Zqm1/WIiRQAAPBGEOrierbQNea/6Wrjv55jhE7WN04xi4niVw0AgC/+OnZxPR3+XWPunBN0jJBXEGqcbRYbxWKKAAD4Igh1ce6KUOXJELrGAowRcgehGIIQAAB+CEJdXPNg6UCzxryvDbSyNF1jAAAEx1/HLq6Ho7GSc7yuOQgF32LDf/o8XWMAAARHEOriYqMbu8aO1zWvLB1O19iJpiDkIAgBAOCHINTFxTXtGH+irvUtNty5yBmgIhQTSRACAMAXQaiLi4327xoLusWGuemq/xgh9/sAAIBmBKEuzj22x6trLMj0+UCbrjbPGuNXDQCAL/46dnHurjF3oJGCD5aOaGn6PF1jAAD4IQh1cc1dYyEMlnZPn286bxiGTjbQNQYAQDAEoS4uLtCssWCDpe3eK0vXOw0zNFERAgDAH0GoizNnjQXsGvO+tnnT1cbnJxuaXxMTza8aAABf/HXs4poHS7e+6arvytLu8UE2mxTtm5oAAABBqKuLNQdLu8xKUNDB0j7T50/WNW2vERlhzigDAADNCEJdXJzHIGd395hZEbK3PH3e3TXGQGkAAAIjCHVxnoOc3QOmXU1jgHy7xnynz7tXo46J5NcMAEAg/IXs4ux2mzlOyD3mp7UtNtxdZ82LKVIRAgAgEIKQBcT5rCXUvI6Q93W+0+fdawgRhAAACIwgZAExPhWhoOsI+UyfN7vG2F4DAICA+AtpAY6mIGN2jYU4a6y2ga4xAABaYpkgVFZWpqlTpyo+Pl6JiYmaOXOmqqurW7z+rrvu0uDBgxUbG6uzzz5bd999tyoqKjqw1W3D0TRgurapqyv4rDF5nXcHp1iCEAAAAVkmCE2dOlWffvqpNmzYoHfeeUebN2/WHXfcEfT6b7/9Vt9++62eeOIJ7dq1Sy+++KLy8/M1c+bMDmx124jxqwg1HvftGgs6a4wgBABAQJGd3YBQ7N69W/n5+froo480evRoSdLy5ct13XXX6YknnlC/fv38XjN06FC98cYb5vNzzz1Xjz76qH7xi1+ooaFBkZGW+OqSJEfT9HezIhSka8w9RsicNdZ0vYMxQgAABGSJv5BFRUVKTEw0Q5AkZWVlyW6368MPPwz5fSoqKhQfH99iCKqtrVVlZaXXo7MFGywd6hYbdI0BABCYJYJQSUmJkpOTvY5FRkaqd+/eKikpCek9jh07pkceeaTF7jRJys3NVUJCgvlIS0s75Xa3Fd+KULDB0r7T50+wjhAAAC3q1CC0YMEC2Wy2Fh979uw57c+prKzU9ddfrwsvvFAPPvhgi9fm5OSooqLCfBw8ePC0P/90uQdL+02fb2X3+dr6xh+oCAEAEFinDpSZP3++br311havGThwoFJTU3XkyBGv4w0NDSorK1NqamqLr6+qqtKECRPUs2dPvfXWW4qKimrxeofDIYfDEVL7O4p7sLTvGCG/LTZ8N12tZx0hAABa0qlBqE+fPurTp0+r12VmZqq8vFzFxcUaNWqUJGnjxo1yuVzKyMgI+rrKykqNHz9eDodDa9euVUxMTJu1vSP5Tp8P1jXmO32erjEAAFpmiVLBBRdcoAkTJmjWrFnaunWr/vGPf2ju3Ln62c9+Zs4YO3TokIYMGaKtW7dKagxB11xzjWpqavT888+rsrJSJSUlKikpkdPp7MyvEzazItTKYGnf6fPsNQYAQMssM4f8lVde0dy5czVu3DjZ7XbddNNNeuqpp8zz9fX12rt3r44fPy5J2r59uzmjbNCgQV7vtX//fqWnp3dY20+Xf0Wo8Xir0+fr2WsMAICWWCYI9e7dW6+++mrQ8+np6ebYGEm68sorvZ5bmd+CisGmzwedNWaJwh8AAB2Ov5AWEHyLDe/r3AUi92DqWtYRAgCgRQQhC/CsCBmGIXehy2+LDXPWWONzBksDANAygpAFeFaE3NUeKdCssaZ1hAzfMUL8mgEACIS/kBbg8KgINXgEoUifFRWDbbFBRQgAgMAIQhYQrCIUaW9l93mCEAAALSIIWUBMkIpQa9Pna5k+DwBAiwhCFuBZEWpwbyQm/8HSntPnnS5DdU72GgMAoCUEIQvwHCPUvM9Yc/Bx85w+7x4fJDFYGgCAYPgLaQExnhWhpiDkO1Ba8p4+7xWEIqkIAQAQCEHIAhwee425K0K+A6Wl5gpRg8ulk02LL0ZH2v0qRwAAoBFByALcg51PNrhU3zTux3egtCRFRbiDkKETdU0zxiL5FQMAEAx/JS3A0RRm6jy7xgIEocimPTcanM1jhGKj6RYDACAYgpAFeE5/r6ltkBRkjFBTOKp3ulTbwBpCAAC0hiBkAQ6P7q3jTV1egSpCUU3hqLFrrGkNIQZKAwAQFEHIAqIi7Ga1p7qpIhRojFCke4yQ09W8vQZdYwAABEUQsgh3VcjsGgtUEbI3V4RONjBYGgCA1vBX0iLcY32qWxgj1FwR8pg1xhghAACCIghZhLuyU91CRcgdhOo91hFiew0AAIIjCFmEu7JT08IYIXfXmGFIx5uuY3sNAACC46+kRTjMIBR81pi7IiQ1V47oGgMAIDiCkEW4KztVJ1sYI2RvPua+jiAEAEBwBCGLcK8H1FLXmGdFiCAEAEDrCEIW4a4I1dS1MFja7tk1Vi+JwdIAALSEIGQRvtPnA1WEbDabGYaaK0L8igEACIa/khbhO2ssUEVIau4eY7A0AACtIwhZhNk15p41FmCwtNQ8hb66qSJE1xgAAMERhCzCEemzsnQrFaGqpuscdI0BABAUfyUtIpQFFaXmSlE1s8YAAGgVQcgi3F1jDS5DUgsVoabjJ5p2n6drDACA4AhCFuHuGnOLsAf+1XmuJSRREQIAoCUEIYvwnQYfHRm4IhTlE5CYPg8AQHD8lbQI38pOdJBZY74VIbrGAAAIjiBkEf4VoSBByK8iRBACACAYgpBFxPiMEQoWhKJ8KkJMnwcAIDj+SlqEf9dY4EqP70KLcdGR7dYmAACsjiBkEb6VneBdY4wRAgAgVJYJQmVlZZo6dari4+OVmJiomTNnqrq6OqTXGoaha6+9VjabTW+//Xb7NrSd+FWEggUhj66x6Eh70IUXAQCAhYLQ1KlT9emnn2rDhg165513tHnzZt1xxx0hvXbZsmWy2awdCEIfI9R8PC6aahAAAC2xxACS3bt3Kz8/Xx999JFGjx4tSVq+fLmuu+46PfHEE+rXr1/Q1+7YsUNLlizRtm3b1Ldv31Y/q7a2VrW1tebzysrK0/8CbcB31pgjyPR5z64wusUAAGiZJSpCRUVFSkxMNEOQJGVlZclut+vDDz8M+rrjx4/r5z//uVasWKHU1NSQPis3N1cJCQnmIy0t7bTb3xZC7RrzvC6WihAAAC2yRBAqKSlRcnKy17HIyEj17t1bJSUlQV83b948XXrppbrxxhtD/qycnBxVVFSYj4MHD55yu9vSqQQhusYAAGhZpwahBQsWyGaztfjYs2fPKb332rVrtXHjRi1btiys1zkcDsXHx3s9ugK/rrEgQYiuMQAAQtepY4Tmz5+vW2+9tcVrBg4cqNTUVB05csTreENDg8rKyoJ2eW3cuFH79u1TYmKi1/GbbrpJY8eOVWFh4Wm0vOOFOljaMzDFsoYQAAAt6tS/lH369FGfPn1avS4zM1Pl5eUqLi7WqFGjJDUGHZfLpYyMjICvWbBggW6//XavYxdffLGefPJJTZw48fQb38HsdpvOiI5QTZ1TUvCtMzyrQHFUhAAAaJElSgYXXHCBJkyYoFmzZikvL0/19fWaO3eufvazn5kzxg4dOqRx48bp5Zdf1pgxY5SamhqwWnT22WdrwIABHf0V2kR8bJQZhM4IUu3xHCDNGCEAAFpmicHSkvTKK69oyJAhGjdunK677jpdfvnlevbZZ83z9fX12rt3r44fP96JrWxf8TFR5s89YwIHIYdHFSg+NirgNQAAoJElKkKS1Lt3b7366qtBz6enp8swjBbfo7XzXd0ZjgiPn4NUhDyCUGIcQQgAgJZYpiIEeW2X4RmKPMV7VIoSqAgBANAigpBFOSIDB6Hk+Bjz515x0R3VHAAALIkgZCEX9m19TaOUeIf581m9YtuzOQAAWJ5lxghBmnn5QG3+4phuvTQ96DWp8TEaf1GKyo/Xa2j/hI5rHAAAFmQzrD6CuJ1VVlYqISFBFRUVXWaVaQAA0LJQ/37TNQYAALotghAAAOi2CEIAAKDbIggBAIBuiyAEAAC6LYIQAADotghCAACg2yIIAQCAbosgBAAAui2CEAAA6LYIQgAAoNsiCAEAgG6LIAQAALotghAAAOi2Iju7AV2dYRiSpMrKyk5uCQAACJX777b773gwBKFWVFVVSZLS0tI6uSUAACBcVVVVSkhICHreZrQWlbo5l8ulb7/9Vj179pTNZmuz962srFRaWpoOHjyo+Pj4Nntf+ONedwzuc8fgPncM7nPHaM/7bBiGqqqq1K9fP9ntwUcCURFqhd1u11lnndVu7x8fH8//yDoI97pjcJ87Bve5Y3CfO0Z73eeWKkFuDJYGAADdFkEIAAB0WwShTuJwOLRo0SI5HI7Obsr3Hve6Y3CfOwb3uWNwnztGV7jPDJYGAADdFhUhAADQbRGEAABAt0UQAgAA3RZBCAAAdFsEoU6yYsUKpaenKyYmRhkZGdq6dWtnN6lL27x5syZOnKh+/frJZrPp7bff9jpvGIYWLlyovn37KjY2VllZWfriiy+8rikrK9PUqVMVHx+vxMREzZw5U9XV1V7X7Ny5U2PHjlVMTIzS0tL0+9//vr2/WpeRm5ur//f//p969uyp5ORkTZo0SXv37vW65uTJk5ozZ47OPPNM9ejRQzfddJNKS0u9rjlw4ICuv/56xcXFKTk5Wb/+9a/V0NDgdU1hYaFGjhwph8OhQYMG6cUXX2zvr9dlPPPMMxo2bJi5gFxmZqbeffdd8zz3uH0sXrxYNptN9957r3mMe902HnzwQdlsNq/HkCFDzPNd/j4b6HCrV682oqOjjVWrVhmffvqpMWvWLCMxMdEoLS3t7KZ1WevXrzd+97vfGW+++aYhyXjrrbe8zi9evNhISEgw3n77beOf//ynccMNNxgDBgwwTpw4YV4zYcIEY/jw4cYHH3xg/P3vfzcGDRpkTJkyxTxfUVFhpKSkGFOnTjV27dplvPbaa0ZsbKzx3//93x31NTvV+PHjjRdeeMHYtWuXsWPHDuO6664zzj77bKO6utq85s477zTS0tKMgoICY9u2bcYPfvAD49JLLzXPNzQ0GEOHDjWysrKMjz/+2Fi/fr2RlJRk5OTkmNd89dVXRlxcnJGdnW189tlnxvLly42IiAgjPz+/Q79vZ1m7dq2xbt064/PPPzf27t1r/Pa3vzWioqKMXbt2GYbBPW4PW7duNdLT041hw4YZ99xzj3mce902Fi1aZFx00UXG4cOHzcfRo0fN8139PhOEOsGYMWOMOXPmmM+dTqfRr18/Izc3txNbZR2+QcjlchmpqanGH/7wB/NYeXm54XA4jNdee80wDMP47LPPDEnGRx99ZF7z7rvvGjabzTh06JBhGIbx9NNPG7169TJqa2vNa+677z5j8ODB7fyNuqYjR44Ykoz33nvPMIzGexoVFWX85S9/Ma/ZvXu3IckoKioyDKMxsNrtdqOkpMS85plnnjHi4+PN+/qb3/zGuOiii7w+a/Lkycb48ePb+yt1Wb169TKee+457nE7qKqqMs477zxjw4YNxhVXXGEGIe5121m0aJExfPjwgOescJ/pGutgdXV1Ki4uVlZWlnnMbrcrKytLRUVFndgy69q/f79KSkq87mlCQoIyMjLMe1pUVKTExESNHj3avCYrK0t2u10ffvihec0Pf/hDRUdHm9eMHz9ee/fu1XfffddB36brqKiokCT17t1bklRcXKz6+nqv+zxkyBCdffbZXvf54osvVkpKinnN+PHjVVlZqU8//dS8xvM93Nd0x//+nU6nVq9erZqaGmVmZnKP28GcOXN0/fXX+90P7nXb+uKLL9SvXz8NHDhQU6dO1YEDByRZ4z4ThDrYsWPH5HQ6vX7hkpSSkqKSkpJOapW1ue9bS/e0pKREycnJXucjIyPVu3dvr2sCvYfnZ3QXLpdL9957ry677DINHTpUUuM9iI6OVmJiote1vve5tXsY7JrKykqdOHGiPb5Ol/PJJ5+oR48ecjgcuvPOO/XWW2/pwgsv5B63sdWrV2v79u3Kzc31O8e9bjsZGRl68cUXlZ+fr2eeeUb79+/X2LFjVVVVZYn7zO7zAPzMmTNHu3bt0pYtWzq7Kd9LgwcP1o4dO1RRUaHXX39d06dP13vvvdfZzfpeOXjwoO655x5t2LBBMTExnd2c77Vrr73W/HnYsGHKyMjQOeecoz//+c+KjY3txJaFhopQB0tKSlJERITfiPnS0lKlpqZ2UquszX3fWrqnqampOnLkiNf5hoYGlZWVeV0T6D08P6M7mDt3rt555x1t2rRJZ511lnk8NTVVdXV1Ki8v97re9z63dg+DXRMfH2+J/9NsC9HR0Ro0aJBGjRql3NxcDR8+XH/84x+5x22ouLhYR44c0ciRIxUZGanIyEi99957euqppxQZGamUlBTudTtJTEzU+eefry+//NIS/00ThDpYdHS0Ro0apYKCAvOYy+VSQUGBMjMzO7Fl1jVgwAClpqZ63dPKykp9+OGH5j3NzMxUeXm5iouLzWs2btwol8uljIwM85rNmzervr7evGbDhg0aPHiwevXq1UHfpvMYhqG5c+fqrbfe0saNGzVgwACv86NGjVJUVJTXfd67d68OHDjgdZ8/+eQTr9C5YcMGxcfH68ILLzSv8XwP9zXd+b9/l8ul2tpa7nEbGjdunD755BPt2LHDfIwePVpTp041f+Zet4/q6mrt27dPffv2tcZ/06c93BphW716teFwOIwXX3zR+Oyzz4w77rjDSExM9BoxD29VVVXGxx9/bHz88ceGJGPp0qXGxx9/bHzzzTeGYTROn09MTDT+93//19i5c6dx4403Bpw+f8kllxgffvihsWXLFuO8887zmj5fXl5upKSkGNOmTTN27dplrF692oiLi+s20+dnz55tJCQkGIWFhV7TYI8fP25ec+eddxpnn322sXHjRmPbtm1GZmamkZmZaZ53T4O95pprjB07dhj5+flGnz59Ak6D/fWvf23s3r3bWLFiRbeabrxgwQLjvffeM/bv32/s3LnTWLBggWGz2Yy//e1vhmFwj9uT56wxw+Bet5X58+cbhYWFxv79+41//OMfRlZWlpGUlGQcOXLEMIyuf58JQp1k+fLlxtlnn21ER0cbY8aMMT744IPOblKXtmnTJkOS32P69OmGYTROoX/ggQeMlJQUw+FwGOPGjTP27t3r9R7//ve/jSlTphg9evQw4uPjjRkzZhhVVVVe1/zzn/80Lr/8csPhcBj9+/c3Fi9e3FFfsdMFur+SjBdeeMG85sSJE8Z//ud/Gr169TLi4uKMH//4x8bhw4e93ufrr782rr32WiM2NtZISkoy5s+fb9TX13tds2nTJmPEiBFGdHS0MXDgQK/P+L677bbbjHPOOceIjo42+vTpY4wbN84MQYbBPW5PvkGIe902Jk+ebPTt29eIjo42+vfvb0yePNn48ssvzfNd/T7bDMMwTr+uBAAAYD2MEQIAAN0WQQgAAHRbBCEAANBtEYQAAEC3RRACAADdFkEIAAB0WwQhAADQbRGEAABAt0UQAtDl3HrrrZo0aVKnff60adP02GOPdchnLViwQHfddVeHfBYAf6wsDaBD2Wy2Fs8vWrRI8+bNk2EYSkxM7JhGefjnP/+pq666St9884169OjR7p937NgxDRw4UDt27NDAgQPb/fMAeCMIAehQJSUl5s9r1qzRwoULtXfvXvNYjx49OiSABHP77bcrMjJSeXl5HfaZN998s9LT0/WHP/yhwz4TQCO6xgB0qNTUVPORkJAgm83mdaxHjx5+XWNXXnml7rrrLt17773q1auXUlJStHLlStXU1GjGjBnq2bOnBg0apHfffdfrs3bt2qVrr71WPXr0UEpKiqZNm6Zjx44FbZvT6dTrr7+uiRMneh1/+umndd555ykmJkYpKSn6yU9+Yp5zuVzKzc3VgAEDFBsbq+HDh+v111/3ev2nn36q//iP/1B8fLx69uypsWPHat++feb5iRMnavXq1adyOwGcJoIQAEt46aWXlJSUpK1bt+quu+7S7NmzdfPNN+vSSy/V9u3bdc0112jatGk6fvy4JKm8vFxXXXWVLrnkEm3btk35+fkqLS3VT3/606CfsXPnTlVUVGj06NHmsW3btunuu+/Www8/rL179yo/P18//OEPzfO5ubl6+eWXlZeXp08//VTz5s3TL37xC7333nuSpEOHDumHP/yhHA6HNm7cqOLiYt12221qaGgw32PMmDH617/+pa+//rqN7xqAVrXJHvYAcApeeOEFIyEhwe/49OnTjRtvvNF8fsUVVxiXX365+byhocE444wzjGnTppnHDh8+bEgyioqKDMMwjEceecS45pprvN734MGDhiRj7969Advz1ltvGREREYbL5TKPvfHGG0Z8fLxRWVnpd/3JkyeNuLg44/333/c6PnPmTGPKlCmGYRhGTk6OMWDAAKOuri7IXTCMiooKQ5JRWFgY9BoA7SOyk3MYAIRk2LBh5s8RERE688wzdfHFF5vHUlJSJElHjhyR1DjoedOmTQHHG+3bt0/nn3++3/ETJ07I4XB4Dei++uqrdc4552jgwIGaMGGCJkyYoB//+MeKi4vTl19+qePHj+vqq6/2ep+6ujpdcsklkqQdO3Zo7NixioqKCvrdYmNjJcmsZgHoOAQhAJbgGyRsNpvXMXd4cblckqTq6mpNnDhRjz/+uN979e3bN+BnJCUl6fjx46qrq1N0dLQkqWfPntq+fbsKCwv1t7/9TQsXLtSDDz6ojz76SNXV1ZKkdevWqX///l7v5XA4JDWHnJaUlZVJkvr06dPqtQDaFkEIwPfSyJEj9cYbbyg9PV2RkaH9X92IESMkSZ999pn5syRFRkYqKytLWVlZWrRokRITE7Vx40ZdffXVcjgcOnDggK644oqA7zls2DC99NJLqq+vD1oV2rVrl6KionTRRReF9R0BnD4GSwP4XpozZ47Kyso0ZcoUffTRR9q3b5/++te/asaMGXI6nQFf06dPH40cOVJbtmwxj73zzjt66qmntGPHDn3zzTd6+eWX5XK5NHjwYPXs2VO/+tWvNG/ePL300kvat2+ftm/fruXLl+ull16SJM2dO1eVlZX62c9+pm3btumLL77Q//zP/3gtGfD3v/9dY8eODal6BKBtEYQAfC/169dP//jHP+R0OnXNNdfo4osv1r333qvExETZ7cH/r+/222/XK6+8Yj5PTEzUm2++qauuukoXXHCB8vLy9Nprr5nVm0ceeUQPPPCAcnNzdcEFF2jChAlat26dBgwYIEk688wztXHjRlVXV+uKK67QqFGjtHLlSq/q0OrVqzVr1qx2uhMAWsKCigDg4cSJExo8eLDWrFmjzMzMdv+8d999V/Pnz9fOnTtD7sID0HaoCAGAh9jYWL388sstLrzYlmpqavTCCy8QgoBOQkUIAAB0W1SEAABAt0UQAgAA3RZBCAAAdFsEIQAA0G0RhAAAQLdFEAIAAN0WQQgAAHRbBCEAANBtEYQAAEC39f8DNqvmcqf8jZIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb0AAAG2CAYAAADr6ViHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/av/WaAAAACXBIWXMAAA9hAAAPYQGoP6dpAABPJ0lEQVR4nO3deXhU9aH/8c/MZLLPJASSCQkJCVsChD2AcW0VWcS1VKvlUopWby2gSGsFf4rXFbW91mpVqr1V6INitWrdQCkoyhXZwxYIWyBhmcmeSUgyk8x8f39wz2kmmSTnTGafz+t5eJ7OmZMzZ04hb79n1QghBIiIiCKANtArQERE5C+MHhERRQxGj4iIIgajR0REEYPRIyKiiMHoERFRxGD0iIgoYjB6REQUMRg9IiKKGIweERFFjLCJ3iuvvIKcnBzExsZi6tSp2LFjR6BXiYiIgkxYRO/dd9/F0qVL8dhjj2HPnj0YN24cZsyYgcrKykCvGhERBRFNONxweurUqZg8eTL+9Kc/AQCcTieysrKwePFiLFu2LMBrR0REwSIq0CvQV3a7Hbt378by5cvlaVqtFtOmTcO2bdvc/ozNZoPNZpNfO51O1NbWon///tBoND5fZyIi8h4hBBobG5GRkQGttucdmCEfverqajgcDphMJpfpJpMJR44ccfszK1euxOOPP+6P1SMiIj+pqKjAoEGDepwn5KPnieXLl2Pp0qXy64aGBmRnZ+O2226DXq8P4JoR+Y8QAkIIOBwOtLe3o6WlBU1NTbDZbNBqtUhPT8fo0aMxceJEZGVloaamBlu2bMHGjRtx8uRJGAwGDB8+HOnp6YiOjkZ7ezscDgf3lpDftbW14e9//zsMBkOv84Z89AYMGACdTgeLxeIy3WKxID093e3PxMTEICYmpst0vV6P6Ohon6wnUTDpHDybzYbGxkZYrVbY7XYkJydjwIAByM/PR2FhIUwmE06fPo0jR47AYDDA4XCgqqoKsbGxiI6ORnp6OuLi4uBwOBg+Chglf+9C/uzN6OhoTJo0CZs2bZKnOZ1ObNq0CUVFRQFcM6Lg1DF4NpsNTU1NqK+vR11dHZqamqDT6WAymTB06FBkZ2dDq9WiqqoK1dXViIuLQ15eHvLz8xEfH4/z58/j2LFjMJvNaGtrg06ng06nQxicH0dhKuRHegCwdOlSzJ8/H4WFhZgyZQpefPFFXLhwAQsWLAj0qhEFFXfBs1qtqKurQ3NzM2JjY5GdnY2CggKMGDECMTExOHXqFOrq6tDQ0AC9Xo/8/Hzo9XpERUWhtLQUZrNZXn56erp8iIAjPgpGYRG9n/zkJ6iqqsKKFStgNpsxfvx4bNiwocvJLUSRTEnwcnJyMH78eBQUFMBoNMJsNuPs2bOora1FdHQ0+vXrh/79+8ujOafTKY/0JNIxPoDho+ATFtEDgEWLFmHRokWBXg2ioKQmeGPGjEG/fv1QW1uLsrIymM1mOBwOGI1G6PV66HQ6JCUlIS8vT14+w0ehImyiR0TuqQne2LFjkZKSgvr6epw4cQJnzpxBS0sLoqOj0dzcLI/wkpKSXMKn1Wrd7upk+CjYMHpEYUxp8CZMmIAxY8a4BK+iogLNzc3QarUQQqC1tVVepkajgdFolMMnBY3ho2DH6BGFKU+CV1dXh5MnT6KiogIXLlyAVquFRqOB0+lEW1tbl+VLI74RI0bI0zuGTwiBgQMHMnwUNBg9ojDkafCOHz+Os2fPysGTbunkdDoBXLwOShrxSXoLHwCGj4IGo0cUZrwdvI6cTifa29sVh+/8+fPya4aPggGjRxRG+ho86RhedzftlUZ8DB+FKkaPKEy4C15DQwPq6+sVB0+j0fR6l3q14Tt69CjDR0GD0SMKA92N8LwdPIma8Gk0Ghw5coTho6DA6BGFODUjvLFjx6Jfv359Cp7E012d0n05GT4KBEaPKIR5Erza2lqcOHGiT8GTeBI+ntVJgcToEYWo7oJXW1uL1tZWxMbGIjc3V77TireDJ1Eavu4uYGf4yJ8YPaIQFCzBkygJ3/Dhw+VpnS9gz8jIkJ9xyfCRLzF6RCEm2IIn6Uv4ADB85BeMHlEIURu85ORkvwRP4mn4pJNbGD7yNUaPKET0FLyWlhbExcUhNzcXEydOREFBgd9GeJ2pCZ90OYPFYpHnYfjIlxg9ohCgNnjJycmoqanByZMn/Ro8SW/hMxqNLiM+ho/8hdEjCnKeBK+2tjZgwZMwfBSMGD2iIKYmeGPGjEFSUlKX4EmPBwoEb4VPo9Ggvb0dABg/6hNGjyhIhXrwJD2Fz2g0qhrxtbe3yw+xJfIEo0cUhPoavJaWlqAInkRt+EpLSxk+8glGjyjIeBq8EydO4Ny5c2hpaYFGowm6KDB8FAwYPaIgojR4kyZNQkFBQcgET6I0fJ1vWeZ0OiGEQGZmJsNHfcLoEQUJT4JXXV2NsrKykAieREn4hg0bJk8rLS1FZWWl/DozMxOxsbHyMhg+UoPRIwoCkRI8CcNHgcLoEQVYpAVPwvBRIDB6RAEUqcGTMHzkb4weUYBEevAkDB/5E6NHFABKg1dYWIjRo0eHbfAk3YVPCCHfuYXhI29g9Ij8TE3wCgoKYDQawzp4ku7Cp9Fouoz4tFotDh8+zPCRaowekR95GryTJ0/i/PnzYRs8idoRX8fwCSEwaNAgho96xOgR+UlPwWtubkZ8fHxEB0+i5CbV3YUPAMNHPWL0iPxASfCGDh0qPx4oUoMnURo+aXscOXKE4SNFGD0iH1MavEmTJmH06NERHzyJkvANHTpUnsbwkRKMHpGPqQmewWBg8DrwNHxCCAAMH3XF6BH5kNPpVBW8mpoaOXitra0RHTyJJ+GrqqqSXzN81BGjR+QjUvDsdrviEV5ZWZkcPIBPCZeoCZ9Go8Hhw4fdhk+j0aCtrY3hi2CMHpEPdAxeY2Mjg+cFvd25xWAwuIz4ugsfAIYvgjF6RF7G4PlOX8InhEBWVhbDF+EYPSIv6hy8uro61NfXdwleQUEBEhMTGTwP9HXEx/BFNkaPyEt6C96wYcPkEZ4UvJMnT8JsNjN4KnkaPumsToYvcjF6RF7A4Pmf2vAdOXIE1dXV8muGLzIxekR95C54dXV1aGlp6XIMj8HzLqXhk7bv4cOHGb4Ix+gR9UFvwes8wquqqkJZWRmD50VKwjdkyBAA/76cgeGLXIwekYcYvOChJnwAR3yRjNEj8gCDF3wYPlKC0SNSqbvgNTc3IyEhgcELIE/C1/E6vuzsbMTHxwNg+MIVo0ekgpLgFRYWYtSoUQxegPR1xMfwhTdGj0ghBi90eDrik67jY/jCF6NHpIDa4FVWVuLUqVMMXgBxxEfuMHpEvVATvISEBAYviKi9nKGkpIThC3OMHlEP1AavqqpKDp7NZgPA4AWa2hFfd+HTaDSw2+0MX4hj9Ii60Vvwhg8fLp+lGR8fLx/Ds1gsDF6Q8eauTrvdDqfTCa1W64c1J29j9IjcUBI8aYTXXfAouPQlfE6nEzk5OQxfGGD0iDph8MJXX0d8DF/oY/SIOmDwwl9P4TMYDHL4Ot+kWrqcgeELbYwe0f9h8CJHd+ETQsgjvtzcXHn64cOHUVNTI79m+EIXo0cE98Grra1FS0uL2+BJlyVUVlYyeCFKya5Ohi/8MHoU8foSvM6/LCm09Ba+xMTEXsOXkJAAgOELFYweRTQlwZs8eTJGjhzJ4IUpb4z4GL7QwehRxOoteCNGjEBhYSGDFwE44oscjB5FJAaPOmP4IgOjRxGnu+A1NzcjMTHRbfDKyspQVVXF4IU5teE7cuQIampqXC5nYPiCG6NHEUVp8EaNGoWYmBhYLBacPn2awYsgnoSvtrYWR48eBfDv8Gk0GthsNoYvyAT0/4mVK1di8uTJMBgMSEtLw80334zS0lKXeVpbW7Fw4UL0798fiYmJmDNnDiwWi8s85eXlmD17NuLj45GWloYHH3wQ7e3t/vwqFALUjPBiYmLkm0dzl2bkcTqdcDqdcviamprQ0NAAq9WKtrY2OXzjx49Hfn4+kpOT5fCdOnUKFy5cgF6vR0xMDDQajRxSCryARm/Lli1YuHAhvv/+e2zcuBFtbW2YPn06Lly4IM/zwAMP4JNPPsF7772HLVu24Ny5c/jRj34kv+9wODB79mzY7XZ89913WL16Nd566y2sWLEiEF+JgpSS4E2ePBmjRo1CbGysHDyO8CJX5/BduHCB4QsDGiHtjA4CVVVVSEtLw5YtW3DllVeioaEBqampePvtt/HjH/8YwMVdCSNHjsS2bdtwySWXYP369bj++utx7tw5mEwmAMCqVavw0EMPoaqqCtHR0b1+rtVqRVJSEubOnatofgotSoPXeYTH4BEAaLVaaLVaREVFITY2FomJiTAajTAajdDr9WhqakJZWRn279+PkpIS1NfXIyUlBSNGjJB3dba1tcFms0EIwV2dPmC327F27Vo0NDTAaDT2OG9Qbf2GhgYAQEpKCgBg9+7daGtrw7Rp0+R58vPzkZ2djW3btgEAtm3bhjFjxsjBA4AZM2bAarXi0KFDbj/HZrPBarW6/KHwxOBRXynd1Tl27FiMGjXKZcRXVlbGEV+QCZoTWZxOJ5YsWYLLLrsMBQUFAACz2Yzo6GgkJye7zGsymWA2m+V5OgZPel96z52VK1fi8ccf9/I3oGDD4JG3qD25paSkBLW1tTh27BgAIDc3Vz6rkye3BFbQRG/hwoU4ePAgtm7d6vPPWr58OZYuXSq/tlqtyMrK8vnnkv/0Fry8vDyXk1YqKyt5lib1yNPwSWd1MnzBISiit2jRInz66af45ptvMGjQIHl6eno67HY76uvrXUZ7FosF6enp8jw7duxwWZ50dqc0T2cxMTGIiYnx8regYOEueNXV1WhtbWXwqE88CV9dXZ1L+BITEwEwfIES0K0thMCiRYvw4YcfYvPmzS5/WQBg0qRJ0Ov12LRpkzyttLQU5eXlKCoqAgAUFRXhwIEDqKyslOfZuHEjjEYjRo0a5Z8vQkGDwSNfU3pWZ8djfFL4ysrK0NTUhOjoaB7jC5CAjvQWLlyIt99+G//85z9hMBjkY3BJSUmIi4tDUlIS7rrrLixduhQpKSkwGo1YvHgxioqKcMkllwAApk+fjlGjRmHevHl4/vnnYTab8cgjj2DhwoUczUUYBo/8hSO+0BXQ6L322msAgB/84Acu09988038/Oc/BwD84Q9/gFarxZw5c2Cz2TBjxgy8+uqr8rw6nQ6ffvop7r33XhQVFSEhIQHz58/HE0884a+vQUGAwSN/UxM+jUbD8AWJgEZPySWCsbGxeOWVV/DKK690O8/gwYPx+eefe3PVKIQweBQoPYXPYDD0OOITQmDIkCEMn58FxYksRJ5i8CjQPA2fdDmDFD6NRoPW1laGz8cYPQpZ3QWvpaUFBoOBwSO/6cuID4DLiI/h8y1Gj0JSb8HLz8/H5MmTkZeXJz8toby8nMEjn+lr+HJzc2EwGABcDJ/D4YBOp/PT2kcORo9CjifBO336tLzbk8hXGL7gx+hRSBFCKApefn4+oqOjGTzyO2/v6mT4vIvRo5AhhEB7ezuDR0FPafik43Yc8fkPo0chgSM8CjW9hS8hIQGDBw8G8O/r+KR7dUqXMzB83sfoUdCTgmez2Rg8CilKwyfNd+jQoS6XMzB83sXoUVBTErwpU6YgLy+PwaOgpCR8nXd1dnws0ZAhQ+QHozJ8fcfoUdDqLXgjR46Uz9LU6/Uwm80oLy9n8CjoqNnVCTB8vsToUVBSGzzpOjwGj4KVp+HreFYnw9d3jB4FHQaPwpW3wqfRaNDS0sLweYDRo6DC4FG46y58QggYjcYew+d0OjFs2DB5xMfwqcfoUdBg8ChS9PZYou7Cd/z4cQBg+PqA0aOgwOBRpFETPo1Gg0OHDjF8XsDoUcB1Dl5tbS1qampcgjdlyhSMGDFCPkuzoqKCwaOQp3bEx/D1HaNHAcXgUaRj+PyL0aOA6S54zc3NMBqNDB5FDE+O8dXU1DB8HmD0KCAYPCJXnoZPuoCd4VOG0SO/UxO8qKgoBo8ihifX8TF86jB65Fdqg2exWOTg2e32QK8+kc9xxOdbjB75jbvgSZclSMGbOnUqhg8f7jZ40i8DonCndMSn0WgAuA9fUlISANfwabXaiP93xOiRXzB4ROooCV92djaAf1/H11v4AER8+Bg98rnegjd69GhMnjxZDp7ZbMaZM2cYPIp4asIHoEv4hg4diuTkZAAMn4TRI5/yJHgVFRWoqalh8IjA8Hkbo0c+w+AReYe3w9fe3o6oqKiIDB+jRz7B4BF5l9rwlZSUoLq6GkePHoUQAsOGDWP4wOiRDzB4RL7hafg6ntzSOXw6nQ46nS5i/t0xeuRVSoI3ZcoUDBs2jMEj8oC3wycd44uU8DF65DVKgzd8+HBotVoGj8hDvYUvPj6e4esGo0deoTZ40nV4DB6RZxg+zzB61GfdBa+5uRlJSUkoKCiQj+FJI7wzZ84weER9pCZ80gXs7sKn1Wpx4cKFiAgfo0d9oiR40jE8Bo/I+9SO+DqGTwiB4cOHIzk5GQkJCRERPkaPPKZ0hMfgEfmWp+GTnscXSeFj9Mgj7oJXVVWFlpYWjvCIAsCb4WtubkZ7ezuA8Asfo0eqdQ5eTU2NfNIKg0cUOJ6c3FJVVSUf45PCFx8fH7bhY/RIFQaPKLgxfD1j9Eix3oI3ZswY+RieRqNh8IgCRGn4Oj6Pr3P4+vXrF5bhY/RIEbXBs1gsDB5RACkJX1ZWljwtUsLH6FGvPBnhnT17lsEjCjA14ZOu45PC53Q6MWLECKSkpIRV+Bg96pGS4E2ZMgVDhw5l8IiCkNoRnxS+EydOAIAcPulyhlAPH6NH3WLwiMIDw/dvjB65xeARhZeewpeYmBgx4WP0qIuegtevXz/5TisMHlFoURu+kpISVFZWhlX4GD1yoeQ6PAaPKHSpDZ8QIqxGfIweydwFr7KyEjabDUlJSRg7diyDRxQGPB3xSbcsC+XwMXoEoOfgdd6lCYDBIwpx3YVPCOH25BYAYRE+Ro8UBW/KlCkYMmQIAMBisTB4RGFA7VmdQOiHj9GLcN0Fr7W1FSkpKV2CZzabce7cOQaPKExEWvgYvQjWW/CkO610Dl5dXR2DRxRGIil8jF6E6kvwbDZb0P1FJqK+8Vb4EhMT0dTUhLa2NgDBFz5GLwJ5EryzZ8+ivr4eNptNfqoyEYUXX4RPCIGoqKigCR+jF2EYPCLqSV/CJ4TAiBEj0L9//y4jvmAJH6MXQZQEb8qUKcjNzQXA4BFFqt7CFxcXh6ysLPl5fAC63LklWMPH6EUIBo+I1FASvkGDBrlM77yrMxjDx+hFgM7Bq6qqQnV1tdvgCSHk6/AYPKLIpjZ80gOkgzl8jF6Y6y140q3FpOBJZ2kyeEQEqAufVqsFgKAOH6MXxqTg2e12RcGzWCwMHhF14cmuzs7hGzBgQFCEj9ELUx2DZ7VaOcIjoj7pS/iEEACAAQMGwGAwoLGxMWDhY/TCEHdpEpEveBo+6axOIPDhY/TCjNpdmgweEanh6VmdSsLXcfm+4lH0ysvLcfr0aTQ3NyM1NRWjR49GTEyMt9eNVOptl+a4ceMwefJk5OTkMHhE5DFfjvi0Wq1Pw6c4eqdOncJrr72GdevW4cyZM/I+WgCIjo7GFVdcgXvuuQdz5syRz+Ah/2HwiMifQjV8iup03333Ydy4cSgrK8NTTz2FkpISNDQ0wG63w2w24/PPP8fll1+OFStWYOzYsdi5c6fqFXn22Weh0WiwZMkSeVpraysWLlwon+46Z84cWCwWl58rLy/H7NmzER8fj7S0NDz44IPyHb4jRU/B69+/v0vwnE4ng0dEXuF0OuF0OuXwXbhwAY2NjXLEpPCNHj0ao0ePhslkQktLC06cOIHS0lJUVVVBp9PBYDBAr9fD4XDIv799NXhSNNJLSEjAyZMn0b9//y7vpaWl4eqrr8bVV1+Nxx57DBs2bEBFRQUmT56seCV27tyJP//5zxg7dqzL9AceeACfffYZ3nvvPSQlJWHRokX40Y9+hP/93/8FADgcDsyePRvp6en47rvvcP78efzsZz+DXq/HM888o/jzQ1l3wWtpacGAAQPkY3hS8HhZAhF5kzeO8aWmpvptxKcRHfdTBkBTUxMmTpyIV199FU899RTGjx+PF198EQ0NDUhNTcXbb7+NH//4xwCAI0eOYOTIkdi2bRsuueQSrF+/Htdffz3OnTsHk8kEAFi1ahUeeughVFVVITo6WtE6WK1WJCUlYe7cuYp/Jhj0FDx3IzwGj4h8RavVQqvVIioqCrGxsUhISIDBYJBHcS0tLThz5gxKSkpw8OBBVFZWIjY2FkOGDEF+fj5SU1PhcDjk8Ol0OsUnt9jtdqxduxYNDQ0wGo09r6fXvrGHFi5ciNmzZ2PatGku03fv3o22tjaX6fn5+cjOzsa2bdsAANu2bcOYMWPk4AHAjBkzYLVacejQoW4/02azwWq1uvwJNe6CV1VVJY/wGDwi8qeednXa7XZ5xDdq1CgUFBQgLS0Nra2tOHnyJI4cOeK3XZ2qz96sqanBihUr8NVXX6GysrJLgWtraxUva926ddizZ4/bY4BmsxnR0dFITk52mW4ymWA2m+V5OgZPel96rzsrV67E448/rng9g013wWttbXW7S/P8+fMwm80MHhH5VE+7OhMTE7vs6pRGfCdPnpSn+XpXp+rozZs3D8ePH8ddd90Fk8nk8mgJNSoqKnD//fdj48aNiI2N9WgZnlq+fDmWLl0qv7ZarV2eDxWslARvypQpGDx4MINHRH7nrfAZjUZYrVavh0919L799lts3boV48aN69MH7969G5WVlZg4caI8zeFw4JtvvsGf/vQnfPHFF7Db7aivr3cZ7VksFqSnpwMA0tPTsWPHDpflSmd3SvO4ExMTE5LXFTJ4RBQKgjl8qqOXn5+PlpYWjz9Qcs011+DAgQMu0xYsWID8/Hw89NBDyMrKgl6vx6ZNmzBnzhwAQGlpKcrLy1FUVAQAKCoqwtNPP43KykqkpaUBADZu3Aij0YhRo0b1eR2DCYNHRKHEF+ETQkCv1/cpfKqj9+qrr2LZsmVYsWIFCgoKoNfrXd7v7cwZicFgQEFBgcu0hIQE9O/fX55+1113YenSpUhJSYHRaMTixYtRVFSESy65BAAwffp0jBo1CvPmzcPzzz8Ps9mMRx55BAsXLgzJkVx31ATP4XDAbDYzeEQUcL4IX1tbm0v41F6AoDp6ycnJsFqtuPrqq12mCyGg0Wi8+gv2D3/4A7RaLebMmQObzYYZM2bg1Vdfld/X6XT49NNPce+996KoqAgJCQmYP38+nnjiCa+tQ6D1Frzx48ejsLCQwSOioKQ0fNL5ISUlJTCbzW7D19TUBJvN1qfwqY7e3Llzodfr8fbbb/fpRBZ3vv76a5fXsbGxeOWVV/DKK690+zODBw/G559/7rV1CCYMHhGFAyXhy8zMdImXu/AlJiYCQJfwqemQ6ugdPHgQe/fuRV5entofJRUYPCIKJ2pHfIBr+IQQSEtL6zZ8SqmOXmFhISoqKhg9H1ISvMmTJyM7O5vBI6KQoXTEB0COn9lsdrllmbvwqaE6eosXL8b999+PBx98EGPGjOlyIkvn+2eSOu6CV1lZCZvNxuARUchTEz6Ju/AZDAYAF8On5iEDqqP3k5/8BABw5513ytM0Go1PTmSJNAweEUWC3sIXGxvbbfik434dw2e32xV/turolZWVqf0RUkAIAafTyeARUUTwdMTX8eQWKXw+jd7gwYPV/gj1QgqedCPszsGbMGECCgsLkZ2djfb2dlgsFgaPiEKepyO+zuGTjvEpoSh633//vXxBeG+am5tRVlaG0aNHK16JSNZT8FJTU+WzNBk8IgpHUvgcDoei8Gk0Gpw/f94lfElJSYo/T9F5nvPmzcOMGTPw3nvv4cKFC27nKSkpwcMPP4yhQ4di9+7dilcgknUXvNbWVgaPiCKG0+mEw+GQwyc9lqipqQltbW1y+EaNGoVRo0YhPT0dNptNfixRdXW14s9SNNIrKSnBa6+9hkceeQQ//elPMWLECGRkZCA2NhZ1dXU4cuQImpqacMstt+DLL7/EmDFjPP7ykaKn4KWlpcnH8LKyshg8Igp7He+lqXbEp+b3oeonp+/atQtbt27F6dOn5QeWTpgwAT/84Q+RkpKiZlFBw99PTmfwiIjc02q10Ol00Ol0iI2NRXx8vMsT2FtbW3H27FmUlJSgpKQE58+fR3R0NOrr6xU9Od2ji9MLCws9/kKRrrfgSSetMHhEFIl6GvEZDAa3J7ecO3dO8fJVR488x+AREfVOTfg0Gg2cTqfiy+kYPT9xFzyLxQKbzeYSvEGDBjF4RBTx1ISvra1NcfSU36WTPKYmeA6Hg8EjIkLXszqbm5vR2NiIxsZGl7M61dwLmiM9H/NkhGexWBg8IiIoG/ENHDhQ8fL6FL3W1lbExsb2ZRFhrXPwKisr5QvP09LSMHHiREyaNInBIyLqgZLwKaV696bT6cSTTz6JzMxMJCYmylfFP/roo/if//kftYsLW70Fb8KECQweEZFCPe3qtNlsipejOnpPPfUU3nrrLTz//PMu17QVFBTgL3/5i9rFhSUlI7zJkyczeEREKnQXvubmZsXLUB29NWvW4PXXX8fcuXOh0+nk6ePGjcORI0fULi7sKAleYWEhMjMzGTwiIpXcha+pqUnxz6s+pnf27FkMGzbM7YqofYJtuGHwiIh8r/MxPukp60qoHumNGjUK3377bZfp77//PiZMmKB2cWHD0+A1NDTAbrczeEREKnQc8ak5pqd6pLdixQrMnz8fZ8+ehdPpxAcffIDS0lKsWbMGn376qdrFhQV3wbNYLLDb7b0GT+2j7omI6KKOjyVSSvVI76abbsInn3yCf/3rX0hISMCKFStw+PBhfPLJJ7j22mvVLi7k9Ra8SZMmMXhERD4ijfiU8ug6vSuuuAIbN2705EfDipLgTZo0CZmZmWhra5PfZ/CIiLyn4zG+3qge6e3cuRPbt2/vMn379u3YtWuX2sWFrO6CJx3D6zjCk4J3/vx5Bo+IKIBUR2/hwoWoqKjoMv3s2bNYuHChV1Yq2PUUPJPJJAcvIyPDJXjSRZQMHhFRYKiOXklJCSZOnNhl+oQJE1BSUuKVlQpmngTPYrEweEREQUB19GJiYmCxWLpMP3/+PKKiwvv+1b0FT3rArhQ8i8XCXZpEREFEdfSmT5+O5cuXo6GhQZ5WX1+Phx9+OKzP3lQSvEmTJrkEz2w2c4RHRBREVA/Nfv/73+PKK6/E4MGD5YvRi4uLYTKZ8Le//c3rKxgMOgfPYrHIF54PHDhQfloCg0dEFNxURy8zMxP79+/H2rVrsW/fPsTFxWHBggW44447oNfrfbGOQUEIAYfDAbvdDpvNBrvdjqioKBgMBvTr1w+JiYnQarXy++3t7XA4HBBCBHrViYjo/3h0EC4hIQH33HOPt9claGk0Guh0OsTGxiI5OVm+GLK2thZlZWXQ6XTQ6XQoKCiA0WiEyWSC0+mEEAJNTU3ySJGIiALLo+gdO3YMX331FSorK7v8Ml+xYoVXVizYaDQaREVFIS4uDikpKfL0uro6HD9+XH4tha+jpqYmtLW1MXxERAGmOnpvvPEG7r33XgwYMADp6ekud7fWaDRhGz3Affi0Wi1qamoYPiKiEKA6ek899RSefvppPPTQQ75Yn6DX3YjPXfiSkpJcfpbhIyIKLNXRq6urw6233uqLdQkZvYVPOp43ZswYho+IKIiojt6tt96KL7/8Er/85S99sT4ho6fwnTx5Un7N8BERBQ/V0Rs2bBgeffRRfP/99xgzZkyXyxTuu+8+r61csPMkfNIxUIaPiMg7dDqd4nlVR+/1119HYmIitmzZgi1btri8p9FoIip6AEd8RESBpCZ4gAfRKysrU/sjYY/hIyLyP+kaaa1W+R01w/sO0X6kNnwdL/Vg+IiI1JGCp9frXX6f9saj6J05cwYff/wxysvLYbfbXd574YUXPFlkWFASPqfTiXHjxnHER0TkoY7Bi42NVfWEH9XR27RpE2688UYMGTIER44cQUFBAU6dOgUhhNvn7EWa7sJXW1vrMuJj+IiI1OscvMTERFX3fVYdveXLl+M3v/kNHn/8cRgMBvzjH/9AWloa5s6di5kzZ6pdXFhi+IiIvM9d8AwGg2+P6R0+fBjvvPPOxR+OikJLSwsSExPxxBNP4KabbsK9996rdpFhieEjIvKe7oKXkJCA1tZWxctRHb2EhAT5ON7AgQNx4sQJjB49GgBQXV2tdnFhzV34NBpNl5NbGD4iou71FDyNRgOr1ap4Waqjd8kll2Dr1q0YOXIkrrvuOvz617/GgQMH8MEHH+CSSy5Ru7iwp/SszrFjxyI5OdnlZxk+Iop0HYMXFxeHhIQEl+A1NDTg9OnTipenOnovvPACmpqaAACPP/44mpqa8O6772L48OERfeZmTxg+IiL1Oo/w3AWvrKwMhw4dUrxM1dEbMmSI/L8TEhKwatUqtYuISGrC169fP17HR0QRrbddmlLwDhw4gBMnTiherkfR27lzJ/r37+8yvb6+HhMnTnT5BU6u1IYvIyNDnsbwEVGkUBO8Y8eOobGxUfGyVUfv1KlTcDgcXabbbDacPXtW7eIiTk/h63iLN4aPiCKR2uBZrVb069cP58+fV7R8xdH7+OOP5f/9xRdfuJxt6HA4sGnTJuTk5Cj/ZhGstxGfEAIAw0dEkUVN8I4ePYrGxkakpKRgyJAh3o/ezTffDODiL+z58+e7vKfX65GTk4P//u//Vv7tIpxGo4FOp3N7HZ804hNCYNy4cQwfEYU9d2dpGo1GxMfH9xi8/Px8ZGVl4dtvv1X0OYqjJ/2Czc3Nxc6dOzFgwADPvhnJpLsIuAvfqVOn5NcMHxGFs+7O0lQSvKFDhyI6OlrxZ3nl0UL19fVdTrUnZTwJHx9ES0ThwpMRXv/+/ZGXl4ehQ4fCYDCoujhd+Q3L/s9zzz2Hd999V3596623IiUlBZmZmdi3b5/axREuhq/jrs60tDSkpKTA6XSirKwMO3fuxL59+1BXVweDwYCBAwciLS1NvtGqmvvOEREFC28Er7W1Vb52XAnVvy1XrVqFrKwsAMDGjRvxr3/9Cxs2bMCsWbPw4IMPql0c/R934evfvz+EEDh16hTDR0RhRUnwTp48qSh47q4o6I7q3Ztms1mO3qefforbbrsN06dPR05ODqZOnap2cdRBd7s6a2pq3O7q7Ii7OokoVLg7hucueAcPHlQUPJ8+T69fv36oqKhAVlYWNmzYgKeeegrAxTMN1dSW3FMSPiEExo8fz/ARUchROsJTEzyfjvR+9KMf4ac//SmGDx+OmpoazJo1CwCwd+9eDBs2TO3iyI3ewiddxzdu3DiX9wGGj4iCl6fBy8/Px5AhQ9wGT6fToa2tTfE6qI7eH/7wB+Tk5KCiogLPP/88EhMTAQDnz5/Hr371K7WLo270FL6OdxRn+IgoFPgqeE6nUx4IKKE6enq9Hr/5zW+6TH/ggQfULop6oSZ8ne+FyvARUbDoLnjSnVbq6+tdTlppamqSgzd06FAkJib2GLyON+jvjaLoffzxx5g1axb0er3L7cjcufHGGxV/OPWO4SOiUOZJ8AYMGCAfw0tMTERzczOam5v7HDxAYfRuvvlmmM1mpKWlybcjc0ej0ag+meXs2bN46KGHsH79ejQ3N2PYsGF48803UVhYCODiSRuPPfYY3njjDdTX1+Oyyy7Da6+9huHDh8vLqK2txeLFi/HJJ59Aq9Vizpw5+OMf/yjveg11PV3AzvARUbAKtuABCqPX8RemN3951tXV4bLLLsMPf/hDrF+/HqmpqTh27JjLWYnPP/88XnrpJaxevRq5ubl49NFHMWPGDJSUlCA2NhYAMHfuXJw/fx4bN25EW1sbFixYgHvuuQdvv/2219Y10KTwxcbGdhs+6axOho+IAq2vwYuPj/d68AAPjul503PPPYesrCy8+eab8rTc3Fz5fwsh8OKLL+KRRx7BTTfdBABYs2YNTCYTPvroI9x+++04fPgwNmzYgJ07d8qjw5dffhnXXXcdfv/737vcrzLUabVaREVF9Rg+AAwfEQVUd9fhqQlea2ur14MHqLwji9PpxF//+ldcf/31KCgowJgxY3DjjTdizZo1qs6ekXz88ccoLCzErbfeirS0NEyYMAFvvPGG/H5ZWRnMZjOmTZsmT0tKSsLUqVOxbds2AMC2bduQnJwsBw8Apk2bBq1Wi+3bt7v9XJvNBqvV6vInVHQOX2pqKlJSUiCEwOnTp7Fr1y4UFxejpqYGRqMRmZmZvHMLEflN5xFeYmKiouBJJ634MniAiugJIXDjjTfiF7/4Bc6ePYsxY8Zg9OjROH36NH7+85/jlltuUf3hJ0+elI/PffHFF7j33ntx3333YfXq1QAu3v0FAEwmk8vPmUwm+T3pWGNHUVFRSElJkefpbOXKlUhKSpL/SHeYCRW9hW/Hjh1y+BITE13CFxUVxfARkU94uktTuiyhc/D0er1Xgweo2L351ltv4ZtvvsGmTZvwwx/+0OW9zZs34+abb8aaNWvws5/9TPGHO51OFBYW4plnngEATJgwAQcPHsSqVau6PLPPm5YvX46lS5fKr61Wa8iHT1JbW4uKigr5tbSrMzMzU/4L09jYiPb2du7qJCKvURK8EydOyNfhKQmeVqv1avAAFSO9d955Bw8//HCX4AHA1VdfjWXLlmHt2rWqPnzgwIEYNWqUy7SRI0eivLwcAJCeng4AsFgsLvNYLBb5vfT0dFRWVrq8397ejtraWnmezmJiYmA0Gl3+hKKeRnwVFRXYsWMH9u7di+rqaiQmJiIjIwMmkwkGg4EjPiLyGqXB6zjCS01N9XvwABXR279/P2bOnNnt+7NmzVL9aKHLLrsMpaWlLtOOHj2KwYMHA7h4Ukt6ejo2bdokv2+1WrF9+3YUFRUBAIqKilBfX4/du3fL82zevBlOpzMiboDdW/h27tyJ4uJiho+IfEJN8I4dOyYHLy8vz+/BA1Ts3qytre1ybK0jk8mEuro6VR/+wAMP4NJLL8UzzzyD2267DTt27MDrr7+O119/HcDF6/6WLFmCp556CsOHD5cvWcjIyJCvFxw5ciRmzpyJu+++G6tWrUJbWxsWLVqE22+/PazO3OxJ5/BJIXO3q3PAgAEu24W7OonIU94K3oULF+B0On0ePEBF9Hp7fINOp0N7e7uqD588eTI+/PBDLF++HE888QRyc3Px4osvYu7cufI8v/3tb3HhwgXcc889qK+vx+WXX44NGzbI1+gBwNq1a7Fo0SJcc8018sXpL730kqp1CXUdw9fxKfYMHxH5gruzNA0GQ7e7NC9cuBDw4AGARii81kCr1WLWrFmIiYlx+77NZsOGDRtC8vFCVqsVSUlJmDt3LqKjowO9On3idDrR3t6O5uZm1NfXo6qqSh6BZ2VlYfLkyXL4mpqacO7cOVgsFoaPiBTzdISXn5+P3Nxcl+AJIeRDLZ4Gz263Y+3atWhoaOj1HA3FIz0lZ1OqOXOTfEMa8cXHx8uvgYt3v3E34svMzJSnMXxE1Bu1J61IIzxfBU8txdHreNcUCm4dw9fxRJWO4RNCYMKECQwfESnm6S7NYAkeEODbkJHvdHeMr/OIj+EjIiXCIXgAoxfWOoav435uho+I1FAavH379uH48eO4cOEC0tLSkJeX12vwAPgteACjF/Y6H+OTuAtfamoqw0dELnwdPH9j9CIAw0dEngi34AGMXsRQGr7x48cjLS2N4SOKcOEYPIDRiyh9DZ/dbvfbuhJR4HgavPz8fOTk5ARt8ABGL+K4C59Wq0VNTU2XyxnS0tIwaNAgl6czMHxE4U1J8I4fP479+/d3G7yWlhY0NzcHXfAARi8idTfi63zLMil8HPERRYbuHgAbHx8fFsEDGL2IxfARUUeREDyA0YtovYVPui0rw0cU3roLXkJCAgC4BO/YsWNobm6GyWSST1qJi4sLieABjF7E6yl8Z86ckY/nMXxE4akvwcvJyQmp4AGMHkHZrk6n04lJkyYxfERhpK/BC5Vdmh0xegRA+a5OKXyDBg2S52H4iEKPp8GTjuHFxsaGXPAARo866G1Xp2TixIkwmUwMH1GI8kbwWltbQy54AKNHnTB8ROFNbfBaWlrCJngAo0duMHxE4cmT4HW8LCHUgwcwetQNT8LHO7cQBa+egqfRaFBXV9dr8FpaWtDS0gIAIRk8gNGjHvQUvrNnz8qvpfDxrE6i4NRb8Gpra3H8+HEcOHBAUfB0Ol1IBg9g9KgXDB9RaFM6wouE4AGMHinQ09MZGD6i4KUkeMeOHZNvLSYFb+TIkRg8eHDYBQ9g9Eghd+GTdoucPXtWvo6PJ7cQBQeluzSlY3itra1hHzyA0SMVetrVee7cOfk1w0cUWAxe9xg9UkVJ+IQQmDhxItLT0xk+Ij9j8HrG6JFqakZ8AwcOdAmf1WpFW1ub39aVKJKovSyhc/BiYmLCOngAo0ceUhO+ziM+ho/I+zwZ4aWnpyMvL08OXmtra1gHD2D0qA844iMKDgyecowe9QnDRxRYSoO3b98+HD9+PKKDBzB65AVKT26ZNGmSHD7plmUMH5HnvBG8lpYWtLa2yssL5+ABjB55idIRX8fwSRg+IvUYPM8weuQ1DB+Rf3QXvMTERAAX/811vNOKkuABCPvgAYweeZlGo+n2lmWdj/FlZGQwfEQqeRK8gQMHYsSIET0GL1IweuRV0rE6JeEDwPARqSDtgmTwPMfokdcpCV/H3SgMH1HvpOBFR0czeH3A6JFP9BY+s9mM3bt3yzeqZviIuqcmeMeOHYPNZmPwusHokc8oCV9xcbE8H8NH1JXS4BUXF+PkyZNy8PLz85GVlYXo6GgGrwNGj3xK6a5OjviIulITvBMnTrgNXmtrK4PXAaNHPtdd+KRricxmM/bs2SNPz8jIQFZWlvya4aNI5C54SUlJSEhIAKAueBqNBlqtNsDfKDgweuQXPY34egof79xCkciT4GVkZCAvL4/B6wWjR36jJnxCCGRmZnJXJ0UctcGz2+1dgtfS0gKbzcbgucHokV+pHfExfBRJvDHCY/B6xuiR37kLnzStY/ikk1sYPooESoJ39OhR7Nu3jyO8PmD0KCCUXM6wd+9eeRrDR+HMk+ANHDiQwfMAo0cB4y580uiO4aNI4Wnw8vPzMWjQIAZPJUaPAsqTER8vZ6BwweD5H6NHAdfTMb66ujpYLBaGj8JOX4On1+sZPA8wehQUehvxMXwUTrwRvNbWVgbPA4weBQ2l4et4Vqd0AfuZM2cYPgoJPQVPo9GgpqZGDt7x48fR1tbG4HkRo0dBRUn4iouLAVw86WXQoEHyyS0MHwU7Bi/wGD0KOmrCB6BL+BoaGtDe3u639SVSgsELDoweBaW+hA8Aw0dBxdPgjRw5EpmZmQyeFzF6FLQ8CV/Hk1sYPgoG3QVPejxQTU0Njhw5gv379+PEiRNoa2tDZmYm8vLyGDwfYPQoqDF8FMoYvODD6FHQUxo+6axOho+CgTeC19LSArvdzuB5EaNHIUFJ+Pbt2ydPY/gokBi84MXoUcjoKXz19fUu4RNCICsry+VBtAwf+UPn4BkMBhiNRsXBi4qKYvB8iNGjkKImfACQlZXFszrJb9yN8NQGr7W1lcHzIUaPQg7DR8FIq9X2aZcmg+cfjB6FJE/Cx2N85CtKd2nu27cPJ0+eRHt7OzIzM5Gfn4+MjAwGz48YPQpZSsPndDoBMHzkG54GLy8vzyV4bW1t0Gq18t9r8g1Gj0KakvB1xPCRN3Xepelp8Ox2O4PnJ4wehbyO4YuNjXV5j+EjX9FqtdDpdIqCd+LECTgcDo7wgkBAdxw7HA48+uijyM3NRVxcHIYOHYonn3xSvsgYuHjq+YoVKzBw4EDExcVh2rRpOHbsmMtyamtrMXfuXBiNRiQnJ+Ouu+5CU1OTv78OBZBGo4FGo0F0dDRiY2NhMBgwYMAAJCcnQ6fTwWKx4MCBA9i5cyfKy8uh1+vl+CUlJSEqiv/9R8q5C17Hk1aqq6sVB0/6u0v+EdB/6c899xxee+01rF69GqNHj8auXbuwYMECJCUl4b777gMAPP/883jppZewevVq5Obm4tFHH8WMGTNQUlIi/1f93Llzcf78eWzcuBFtbW1YsGAB7rnnHrz99tuB/HrkZ9Ivjujo6C7vuRvxZWdny9fxlZeXc8RHinQXPOkBsDU1NSgtLWXwglRAo/fdd9/hpptuwuzZswEAOTk5eOedd7Bjxw4AF0d5L774Ih555BHcdNNNAIA1a9bAZDLho48+wu23347Dhw9jw4YN2LlzJwoLCwEAL7/8Mq677jr8/ve/R0ZGRmC+HAWEJ+GTLmdg+Kg3SoLHEV5wC+juzUsvvRSbNm3C0aNHAQD79u3D1q1bMWvWLABAWVkZzGYzpk2bJv9MUlISpk6dim3btgEAtm3bhuTkZDl4ADBt2jRotVps377d7efabDZYrVaXPxQ+lO7q3LVrl7yrc9CgQcjOzkZycjJ3dZJbngRv0KBBDF6QCei/7mXLlsFqtSI/Px86nQ4OhwNPP/005s6dCwAwm80AAJPJ5PJzJpNJfs9sNiMtLc3l/aioKKSkpMjzdLZy5Uo8/vjj3v46FER6G/FVVlZi//798rSOIz5pHo74SNJb8Kqrq7vs0hw0aBBGjBiBjIwM6HQ6Bi9IBDR6f//737F27Vq8/fbbGD16NIqLi7FkyRJkZGRg/vz5Pvvc5cuXY+nSpfJrq9XqcjYfhQdPwtfx7wHDR4Dy4O3duxdlZWUuI7yBAwdCp9PBZrMxeEEioNF78MEHsWzZMtx+++0AgDFjxuD06dNYuXIl5s+fj/T0dACAxWLBwIED5Z+zWCwYP348ACA9PR2VlZUuy21vb0dtba38853FxMQgJibGB9+Igk1v4auqqpLD53Q6kZOTw/CRzBvBa21tRXt7O4MXJAJ6TK+5ubnL7XZ0Op18B43c3Fykp6dj06ZN8vtWqxXbt29HUVERAKCoqAj19fXYvXu3PM/mzZvhdDoxdepUP3wLCnY9HePTarVy+Pbs2YNTp04hKipKvpyBx/giF4MXngL6r/mGG27A008/jezsbIwePRp79+7FCy+8gDvvvBPAxV9WS5YswVNPPYXhw4fLlyxkZGTg5ptvBgCMHDkSM2fOxN13341Vq1ahra0NixYtwu23384zN0mmZsQHgCO+CNfTdXhCCAYvhAU0ei+//DIeffRR/OpXv0JlZSUyMjLwn//5n1ixYoU8z29/+1tcuHAB99xzD+rr63H55Zdjw4YNLnfeWLt2LRYtWoRrrrkGWq0Wc+bMwUsvvRSIr0RBzNPwST9XV1cHh8Phn5WlgFESvCNHjqC4uJjBC0Ea0fH2JxHKarUiKSkJc+fOdfsLkcKLEAJCCNjtdrS2tqKxsRHV1dWor6+H0+lEamoqxo4di4kTJyInJwft7e04c+YMysvLGb4wx+CFJrvdjrVr16KhoQFGo7HHeXmwgiJOdyM+rVaL2tpaecQn/fdgTk6Oy+UMdXV18i81Ch8MXmRg9Cgi9bSrUwrfgQMH5Gmdj/ExfOGFwYscjB5FLOkXk16v7/Iewxc51AbP6XQyeCGM0aOIJj2lmuGLTJ4ET3rieXp6OjQaDYMXYhg9ingMX2RS+nig4uJinDx5EkKILsGz2+0MXohh9IjA8EWa3oJXVVXVJXhZWVkYMWIEgxfiGD2i/6M0fB3P6mT4Qk9fgmcymRi8EMfoEXWgJHwHDx6Uf9ExfKGlr8HTarUMXohj9Ig66Sl80p1bOOILPUqCd/jwYfmkFQYvPDF6RG4oCd/BgwcBXLzDS25uLsMXxBg8kjB6RN1QEz7g4lNBsrOz5dcMX3Bg8KgjRo+oB0p3dUoYvuDiSfCys7MxfPhwBi9MMXpEvegtfNXV1W5HfNKzIhm+wGDwyB1Gj0gBT3Z18ibVgaM2eAAYvAjB6BEppPasToYvMNwFLzk5WX7iubvgScfw0tLSGLwwx+gRqaBmV6cQAkOGDGH4/EhJ8EpKSrBv3z4GL0IxekQq9Ra+mpoal12dDJ9/dBc86ebRHUd4p06dAsDgRSJGj8gDDF9wURK8kpKSHoPX2toKp9Mpn4BE4YnRI/JQd+HTaDSoq6tzG76OF7DX1tbC4XAwfH3kSfCkk1YYvMjD6BH1QU8jPnfhk+7cIu0+q62t5YivDxg8UovRI+ojteHrvKuzpqaGIz4PMHjkCUaPyAuUhk+6nMFd+Nrb2/nLVyEpeHq9vtfglZWVQaPRMHgEgNEj8holJ7ccOnRIntb5Oj6GT5neRniVlZUu1+ExeNQRo0fkRWrC5+46PoavZ72N8DoG79SpUwwedcHoEXmZ2hEfw6eMkhFeSUkJ9u7di9OnT3cJnkajYfCI0SPyBYbPu5Qcwzt8+HCPwbPZbAweMXpEvuJp+KSzOBm+izoGLz4+vtcRnlarRVZWFoNHbjF6RD7kLnxS1Dji6x2DR97G6BH5mLvwSZcuMHzdY/DIFxg9Ij9QuquT1/Fd5K3gCSEiZpuRMowekZ8o2dVZUlIC4OJIcOjQoREZPm8Gj3e5oc4YPSI/UjLik8IHIOLC50nwBg8ejKFDhzJ4pAijR+RnDJ97nYNnNBqRlJTkErxDhw6huLiYwSOPMXpEAeBp+ML1cgZ3I7zOwZPupdk5eKmpqQweKcboEQWIJ+HLzMyUX4dL+Lob4RkMBjidzh5HeKmpqdBqtQweKcboEQWQ0vBJZ3WGW/jcBU86hqckeBzhkVqMHlGA9RS+hoYG1NbW4vDhw9BoNBBCYNiwYXL4NBoNqqurQzJ8aoOn0+mQnZ3tEjy73c7gkSqMHlEQUBK+jiO+juEDEHLh6yl47k5a0el0GDx4MIYNG4b+/fszeOQxRo8oSCgJX8c7t7gLX1tbG3Q6nV/W11O9Bc9iseDQoUPYu3cvysvLGTzyKkaPKIh0Fz6NRoP6+nrU19e7DZ/0yz/Yw8fgUaAxekRBpreTW9yFLyMjQ34drOFj8CgYMHpEQaiv4ausrITD4Qia8DF4FCwYPaIg5a3wabXagIZCSfAOHjyI4uJiBo98jtEjCmJqwiddzhBM4fMkeDk5ORg6dCiDRz7B6BEFuVAd8TF4FIwYPaIQoDR8TqcTQODDpzR4e/fuRUVFBYNHfsPoEYUIJeE7fPiwPE0KnxQNf4Wvt+CZzWb5pJWKigpERUXJtxZj8MjXGD2iEOJJ+AYOHCi/9nX41ASvvLwcer3eJXgAGDzyKUaPKMR4Er6Od27xVfi8Eby2tjYGj3yK0SMKQX0Jn0ajgcVi8Wr4GDwKFYweUYhSEz4hBIYPH+4y4vNW+Bg8CiWMHlEIUzvi83b4ugue9ABYs9mMAwcOYN++fS7BGzZsGFJSUgAweORfjB5RiPM0fFJkPA0fg0ehiNEjCgOehK/jdXwWiwXt7e3Q6XSKAqQ0eHv37sWZM2dcgtevXz8ADB4FBqNHFCb8Fb6+Bk+j0TB4FDCMHlEY8XX41AYvOjoa2dnZDB4FDUaPKMyoPatzxIgRisLnSfBycnIwZMgQBo+CBqNHFIbUjvh6Cp/0TD4Gj8IBo0cUpjwJn7uzOjUaDfR6PYNHYYHRIwpjfRnxaTQaVFVVyU9gT0hIgMFgYPAopDF6RGHO0/BptVpERUWhrq4Oer0eiYmJqoIH8LIECj6MHlEEUBu+vLw8DBo0CLGxsTh37hxaWlqQmJjocmux3oLX3t7O4FHQYfSIIoSS8JWWlgIAoqKiMHr0aAwdOhTp6emoqalBc3MzbDYbqqqq5AfASsHLzc1Fbm4ug0dBj9EjiiA9hc9qtaKhoQHHjh1DbGwskpKSkJ6ejsGDB6OxsRHHjh1DcXEx9u7diwMHDuDcuXOIiYnB4MGDGTwKGYweUYRxFz5pWlNTE1paWlBRUYFjx45h8ODByMrKQkJCAs6dO4fKykocPHgQp06dQnx8vDzCS05OBsDgUfBj9IgiUOfwSdfjxcbGyk8ur6qqwtGjR+XLFKqrq9HY2Ij29nYkJiYiMzMTgwcPRr9+/SCEYPAoJDB6uHhXCuDimWZEkUaj0SAqKgoGgwHx8fFwOp3QaDRobm7Gvn37cPr0aeh0OrS0tKC1tRVDhw6Fw+FAQkICYmJi0NraCofDEeivQRFM+t0t/S7viUYomSvMnTx5EkOHDg30ahARUR9UVFRg0KBBPc7DkR4gP9urvLwcSUlJAV6b4GW1WpGVlYWKigoYjcZAr07Q4nbqHbeRMtxOyggh0NjY6HIrve4werh4I10ASEpK4l8sBYxGI7eTAtxOveM2UobbqXdKByxaH68HERFR0GD0iIgoYjB6AGJiYvDYY48hJiYm0KsS1LidlOF26h23kTLcTt7HszeJiChicKRHREQRg9EjIqKIwegREVHEYPSIiChiRHz0XnnlFeTk5CA2NhZTp07Fjh07Ar1KPvXNN9/ghhtuQEZGBjQaDT766COX94UQWLFiBQYOHIi4uDhMmzYNx44dc5mntrYWc+fOhdFoRHJyMu666y40NTW5zLN//35cccUViI2NRVZWFp5//nlffzWvWblyJSZPngyDwYC0tDTcfPPN8nPmJK2trVi4cCH69++PxMREzJkzBxaLxWWe8vJyzJ49G/Hx8UhLS8ODDz6I9vZ2l3m+/vprTJw4ETExMRg2bBjeeustX389r3nttdcwduxY+cLpoqIirF+/Xn6f28i9Z599FhqNBkuWLJGncVv5kYhg69atE9HR0eKvf/2rOHTokLj77rtFcnKysFgsgV41n/n888/F//t//0988MEHAoD48MMPXd5/9tlnRVJSkvjoo4/Evn37xI033ihyc3NFS0uLPM/MmTPFuHHjxPfffy++/fZbMWzYMHHHHXfI7zc0NAiTySTmzp0rDh48KN555x0RFxcn/vznP/vra/bJjBkzxJtvvikOHjwoiouLxXXXXSeys7NFU1OTPM8vf/lLkZWVJTZt2iR27dolLrnkEnHppZfK77e3t4uCggIxbdo0sXfvXvH555+LAQMGiOXLl8vznDx5UsTHx4ulS5eKkpIS8fLLLwudTic2bNjg1+/rqY8//lh89tln4ujRo6K0tFQ8/PDDQq/Xi4MHDwohuI3c2bFjh8jJyRFjx44V999/vzyd28p/Ijp6U6ZMEQsXLpRfOxwOkZGRIVauXBnAtfKfztFzOp0iPT1d/O53v5On1dfXi5iYGPHOO+8IIYQoKSkRAMTOnTvledavXy80Go04e/asEEKIV199VfTr10/YbDZ5noceekjk5eX5+Bv5RmVlpQAgtmzZIoS4uE30er1477335HkOHz4sAIht27YJIS7+x4VWqxVms1me57XXXhNGo1HeLr/97W/F6NGjXT7rJz/5iZgxY4avv5LP9OvXT/zlL3/hNnKjsbFRDB8+XGzcuFFcddVVcvS4rfwrYndv2u127N69G9OmTZOnabVaTJs2Ddu2bQvgmgVOWVkZzGazyzZJSkrC1KlT5W2ybds2JCcno7CwUJ5n2rRp0Gq12L59uzzPlVdeiejoaHmeGTNmoLS0FHV1dX76Nt7T0NAA4N83Jt+9ezfa2tpctlN+fj6ys7NdttOYMWNgMpnkeWbMmAGr1YpDhw7J83RchjRPKP79czgcWLduHS5cuICioiJuIzcWLlyI2bNnd/k+3Fb+FbE3nK6urobD4XD5SwQAJpMJR44cCdBaBZbZbAYAt9tEes9sNiMtLc3l/aioKKSkpLjMk5ub22UZ0nv9+vXzyfr7gtPpxJIlS3DZZZehoKAAwMXvEB0dLT8tXNJ5O7nbjtJ7Pc1jtVrR0tKCuLg4X3wlrzpw4ACKiorQ2tqKxMREfPjhhxg1ahSKi4u5jTpYt24d9uzZg507d3Z5j3+f/Ctio0ekxMKFC3Hw4EFs3bo10KsSlPLy8lBcXIyGhga8//77mD9/PrZs2RLo1QoqFRUVuP/++7Fx40bExsYGenUiXsTu3hwwYAB0Ol2XM6QsFgvS09MDtFaBJX3vnrZJeno6KisrXd5vb29HbW2tyzzultHxM0LBokWL8Omnn+Krr75yeTBleno67HY76uvrXebvvJ162wbdzWM0GkPmv8qjo6MxbNgwTJo0CStXrsS4cePwxz/+kduog927d6OyshITJ05EVFQUoqKisGXLFrz00kuIioqCyWTitvKjiI1edHQ0Jk2ahE2bNsnTnE4nNm3ahKKiogCuWeDk5uYiPT3dZZtYrVZs375d3iZFRUWor6/H7t275Xk2b94Mp9OJqVOnyvN88803aGtrk+fZuHEj8vLyQmLXphACixYtwocffojNmzd32VU7adIk6PV6l+1UWlqK8vJyl+104MABl/9A2LhxI4xGI0aNGiXP03EZ0jyh/PfP6XTCZrNxG3VwzTXX4MCBAyguLpb/FBYWYu7cufL/5rbyo0CfSRNI69atEzExMeKtt94SJSUl4p577hHJyckuZ0iFm8bGRrF3716xd+9eAUC88MILYu/eveL06dNCiIuXLCQnJ4t//vOfYv/+/eKmm25ye8nChAkTxPbt28XWrVvF8OHDXS5ZqK+vFyaTScybN08cPHhQrFu3TsTHx4fMJQv33nuvSEpKEl9//bU4f/68/Ke5uVme55e//KXIzs4WmzdvFrt27RJFRUWiqKhIfl86xXz69OmiuLhYbNiwQaSmpro9xfzBBx8Uhw8fFq+88kpInWK+bNkysWXLFlFWVib2798vli1bJjQajfjyyy+FENxGPel49qYQ3Fb+FNHRE0KIl19+WWRnZ4vo6GgxZcoU8f333wd6lXzqq6++EgC6/Jk/f74Q4uJlC48++qgwmUwiJiZGXHPNNaK0tNRlGTU1NeKOO+4QiYmJwmg0igULFojGxkaXefbt2ycuv/xyERMTIzIzM8Wzzz7rr6/YZ+62DwDx5ptvyvO0tLSIX/3qV6Jfv34iPj5e3HLLLeL8+fMuyzl16pSYNWuWiIuLEwMGDBC//vWvRVtbm8s8X331lRg/fryIjo4WQ4YMcfmMYHfnnXeKwYMHi+joaJGamiquueYaOXhCcBv1pHP0uK38h48WIiKiiBGxx/SIiCjyMHpERBQxGD0iIooYjB4REUUMRo+IiCIGo0dERBGD0SMioojB6BERUcRg9IgC5Oc//zluvvnmgH3+vHnz8Mwzz/jls5YtW4bFixf75bOIesI7shD5gEaj6fH9xx57DA888ACEEF2eo+YP+/btw9VXX43Tp08jMTHR559XXV2NIUOGoLi4GEOGDPH55xF1h9Ej8gHpwZ4A8O6772LFihUoLS2VpyUmJvolNt35xS9+gaioKKxatcpvn3nrrbciJycHv/vd7/z2mUSdcfcmkQ+kp6fLf5KSkqDRaFymJSYmdtm9+YMf/ACLFy/GkiVL0K9fP5hMJrzxxhu4cOECFixYAIPBgGHDhmH9+vUun3Xw4EHMmjULiYmJMJlMmDdvHqqrq7tdN4fDgffffx833HCDy/RXX30Vw4cPR2xsLEwmE3784x/L7zmdTqxcuRK5ubmIi4vDuHHj8P7777v8/KFDh3D99dfDaDTCYDDgiiuuwIkTJ+T3b7jhBqxbt86TzUnkNYweURBZvXo1BgwYgB07dmDx4sW49957ceutt+LSSy/Fnj17MH36dMybNw/Nzc0AgPr6elx99dWYMGECdu3ahQ0bNsBiseC2227r9jP279+PhoYGFBYWytN27dqF++67D0888QRKS0uxYcMGXHnllfL7K1euxJo1a7Bq1SocOnQIDzzwAP7jP/5Dfkr62bNnceWVVyImJgabN2/G7t27ceedd6K9vV1expQpU3DmzBmcOnXKy1uNSIUAPuGBKCK8+eabIikpqcv0+fPni5tuukl+fdVVV4nLL79cft3e3i4SEhLEvHnz5Gnnz58XAMS2bduEEEI8+eSTYvr06S7LraioEAC6PBJK8uGHHwqdTiecTqc87R//+IcwGo3CarV2mb+1tVXEx8eL7777zmX6XXfdJT9Hcfny5SI3N1fY7fZutoIQDQ0NAoD4+uuvu52HyNeiAtxcIupg7Nix8v/W6XTo378/xowZI08zmUwAID9Be9++ffjqq6/cHh88ceIERowY0WV6S0sLYmJiXE62ufbaazF48GAMGTIEM2fOxMyZM3HLLbcgPj4ex48fR3NzM6699lqX5djtdkyYMAEAUFxcjCuuuAJ6vb7b7xYXFwcA8iiVKBAYPaIg0jkaGo3GZZoUKqfTCQBoamrCDTfcgOeee67LsgYOHOj2MwYMGIDm5mbY7XZER0cDAAwGA/bs2YOvv/4aX375JVasWIH/+q//ws6dO9HU1AQA+Oyzz5CZmemyrJiYGAD/DlpPamtrAQCpqam9zkvkK4weUQibOHEi/vGPfyAnJwdRUcr+OY8fPx4AUFJSIv9vAIiKisK0adMwbdo0PPbYY0hOTsbmzZtx7bXXIiYmBuXl5bjqqqvcLnPs2LFYvXo12trauh3tHTx4EHq9HqNHj1b1HYm8iSeyEIWwhQsXora2FnfccQd27tyJEydO4IsvvsCCBQvgcDjc/kxqaiomTpyIrVu3ytM+/fRTvPTSSyguLsbp06exZs0aOJ1O5OXlwWAw4De/+Q0eeOABrF69GidOnMCePXvw8ssvY/Xq1QCARYsWwWq14vbbb8euXbtw7Ngx/O1vf3O5TOPbb7/FFVdcoWhUSOQrjB5RCMvIyMD//u//wuFwYPr06RgzZgyWLFmC5ORkaLXd//P+xS9+gbVr18qvk5OT8cEHH+Dqq6/GyJEjsWrVKrzzzjvyqOzJJ5/Eo48+ipUrV2LkyJGYOXMmPvvsM+Tm5gIA+vfvj82bN6OpqQlXXXUVJk2ahDfeeMNl1Ldu3TrcfffdPtoSRMrw4nSiCNTS0oK8vDy8++67KCoq8vnnrV+/Hr/+9a+xf/9+xbthiXyBIz2iCBQXF4c1a9b0eBG7N124cAFvvvkmg0cBx5EeERFFDI70iIgoYjB6REQUMRg9IiKKGIweERFFDEaPiIgiBqNHREQRg9EjIqKIwegREVHEYPSIiChi/H8bCkFNLwFxtAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\"\"\"\"Solve 1D elastic equations using displacement\"\"\"\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as pl\n",
    "from math import *\n",
    "\n",
    "#-----------------------------\n",
    "# functions\n",
    "#-----------------------------\n",
    "def ricker(fp,tp,M,dt) :\n",
    "   s=np.zeros((M))\n",
    "   for k in range(0,M):\n",
    "       t=k*dt\n",
    "       wp=2.0*3.14159*fp\n",
    "       s[k] =(1-0.5*pow(wp,2)*pow((t-tp),2))*exp(-0.25*pow(wp,2)*pow(t-tp,2))\n",
    "   return(s)\n",
    "\n",
    "#------------------------------\n",
    "# Main script\n",
    "#-------------------------------\n",
    "M=5000 #No of timesteps\n",
    "N=1000 #No of x grdipoints\n",
    "h=5.0  #Size of each gridcell in the x-direction\n",
    "k=0.0005 #Size of timesteps\n",
    "#k=0.005 #This value leads to unstability\n",
    "c=2000.0 #Wave velocity in m/s\n",
    "r=((k*c)/h) \n",
    "spos=N/2\n",
    "s=ricker(15.0,0.1,M,k)\n",
    "\n",
    "pl.plot(s)\n",
    "pl.xlabel(\"Time (sec)\")\n",
    "pl.ylabel(\"Amplitude \")\n",
    "pl.show()\n",
    "\n",
    "u=np.zeros((N,M))\n",
    "\n",
    "for j in range(1,M-1):\n",
    "  for i in range(1,N-1):\n",
    "    u[i,j+1] = 2*u[i,j] - u[i,j-1] + r*r*(u[i+1,j]- 2*u[i,j]+u[i-1,j])\n",
    "    if i == spos:\n",
    "      u[i,j+1] = u[i,j+1]+s[j]*k*k\n",
    "\n",
    "    \n",
    "#Set aspect ratio\n",
    "pl.imshow(u,cmap='gray')\n",
    "ar=1.0\n",
    "ax=pl.gca()\n",
    "asr = 1.0/(ax.get_data_ratio()*ar)\n",
    "pl.Axes.set_aspect(ax,asr)\n",
    "pl.ylabel(\"Distance (m)\")\n",
    "pl.xlabel(\"Time (sec)\")\n",
    "\n",
    "pl.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15081acf",
   "metadata": {},
   "source": []
  }
 ],
 "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}