{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from matplotlib import animation\n", "from matplotlib.animation import PillowWriter\n", "import numba\n", "from numba import jit" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The TRUE wave equation (with damping and stress-strain coupling):\n", "\n", "$$\\frac{\\partial^2 y}{\\partial x^2} - \\frac{1}{c^2}\\frac{\\partial^2 y}{\\partial t^2} - \\gamma \\frac{\\partial y}{\\partial t} -l^2 \\frac{\\partial^4 y}{\\partial x^4} = 0 $$\n", "\n", "$$ y(0, t) = y(L, t) = 0 $$\n", "$$ y(x, 0) = f(x) $$\n", "\n", "$y$ is the amplitude of the string, $x$ is the location of the string, and $t$ is time.\n", "\n", "Parameters:\n", "\n", "* $c$: speed of the wave $[m/s]$\n", "* $\\gamma$: damping constant $[s/m]$\n", "* $l$: characteristic length (stiffness term) [dimensionless]\n", "\n", "Discrete form:\n", "\n", "$$\\frac{y_{j+1}^{m} -2y_j^m + y_{j-1}^{m}}{\\Delta x^2} - \\frac{1}{c^2}\\frac{y_j^{m+1} -2y_j^m + y_j^{m-1}}{\\Delta t^2} - \\gamma \\frac{y_j^{m+1} - y_j^{m-1}}{2 \\Delta t} - l^2 \\frac{y_{j-2}^m -4y_{j-1}^m +4y_{j}^m -4y_{j+1}^m +y_{j+2}^m}{\\Delta x^4} =0 $$ \n", "\n", "Solve for $y_j^{m+1}$ (the amplitude of the string at the next time):\n", "\n", "$$y_j^{m+1} = \\left[\\frac{1}{c^2 \\Delta t^2} + \\frac{\\gamma}{2 \\Delta t} \\right]^{-1} \\left[\\frac{1}{\\Delta x^2} \\left( y_{j+1}^{m} -2y_j^m + y_{j-1}^{m} \\right) -\\frac{1}{c^2 \\Delta t^2} \\left( y_j^{m-1} - 2y_j^m \\right) + \\frac{\\gamma}{2 \\Delta t}y_j^{m-1} -\\frac{l^2}{\\Delta x^4} \\left( y_{j-2}^m -4y_{j-1}^m +6y_{j}^m -4y_{j+1}^m +y_{j+2}^m \\right) \\right] $$\n", "\n", "In order for this to be stable, need $c \\Delta t/\\Delta x < 1$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Solving" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Guitar string length $\\boxed{L = 0.7~\\text{m}}$\n", "\n", "Choose $\\boxed{N_x=101}$ guitar string positions $\\implies$ $\\boxed{\\Delta x = 0.7~\\text{mm}}$\n", "\n", "Note that the fundamental frequency of a guitar note is $f = c/2L$. With an \"A note\" at 220Hz get $\\boxed{c = 308~\\text{m/s}} $\n", "\n", "To obey our constraint we thus set $\\boxed{\\Delta t = 5 \\times 10^{-6} s}$. In order to get multiple seconds of a result, choose $\\boxed{N_t = 500000}$\n", "\n", "Two parameters that seemed to give a solution that sounded like a string were $\\boxed{l= 2 \\times 10^{-6}}$ and $\\boxed{\\gamma = 2.6 \\times 10^{-5} s/m}$\n", "\n" ] }, { "cell_type": "code", "execution_count": 110, "metadata": {}, "outputs": [], "source": [ "Nx = 101\n", "Nt = 500000\n", "L =0.7\n", "dx = L/(Nx-1)\n", "f = 220\n", "c = 2*L*f\n", "dt = 5e-6\n", "l=5e-5\n", "gamma=5e-5" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Initial state of the string:" ] }, { "cell_type": "code", "execution_count": 111, "metadata": {}, "outputs": [], "source": [ "ya = np.linspace(0, 0.01, 70)\n", "yb = np.linspace(0.01, 0, 31)\n", "y0 = np.concatenate([ya, yb])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create 2D array of $y(x, t)$" ] }, { "cell_type": "code", "execution_count": 112, "metadata": {}, "outputs": [], "source": [ "sol = np.zeros((Nt, Nx))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Make the solution at $t=0$ and $t=1$ equal to the \"pluck\"" ] }, { "cell_type": "code", "execution_count": 113, "metadata": {}, "outputs": [], "source": [ "sol[0] = y0\n", "sol[1] = y0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Go through the iterative procedure:" ] }, { "cell_type": "code", "execution_count": 114, "metadata": {}, "outputs": [], "source": [ "@numba.jit(\"f8[:,:](f8[:,:], i8, i8, f8, f8, f8, f8)\", nopython=True, nogil=True)\n", "def compute_d(d, times, length, dt, dx, l, gamma):\n", " for t in range(1, times-1):\n", " for i in range(2, length-2):\n", " outer_fact = (1/(c**2 * dt**2) + gamma/(2*dt))**(-1)\n", " p1 = 1/dx**2 * (d[t][i-1] - 2*d[t][i] + d[t][i+1])\n", " p2 = 1/(c**2 * dt**2) * (d[t-1][i] - 2*d[t][i])\n", " p3 = gamma/(2*dt) * d[t-1][i]\n", " p4 = l**2 / dx**4 * (d[t][i+2] - 4*d[t][i+1] + 6*d[t][i] - 4*d[t][i-1] + d[t][i-2])\n", " d[t+1][i] = outer_fact * (p1 - p2 + p3 - p4)\n", " return d" ] }, { "cell_type": "code", "execution_count": 115, "metadata": {}, "outputs": [], "source": [ "sol = compute_d(sol, Nt, Nx, dt, dx, l, gamma)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot the string at some sample points" ] }, { "cell_type": "code", "execution_count": 116, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAABLcAAADjCAYAAAB3qPb8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABfKElEQVR4nO3dd3hsV332/e9SLzOj3svpvdocN2xsgynGQAyht1DjkCch5U0BkicE8ry8IaTBE0hx6CUBAqGEDgbbgMG4nd67eteot5nf+8fakkaypCPpSEft/lzXvvaePWvPXnPsOUdza63fcmaGiIiIiIiIiIjISpS01B0QERERERERERGZL4VbIiIiIiIiIiKyYincEhERERERERGRFUvhloiIiIiIiIiIrFgKt0REREREREREZMVSuCUiIiIiIiIiIitWylJ3YLUpLCy09evXL3U3RERERERERERWjSeeeKLVzIqmek7h1gJbv349jz/++FJ3Q0RERERERERk1XDOXZruOU1LFBERERERERGRFUvhloiIiIiIiIiIrFgKt0REREREREREZMVSuCUiIiIiIiIiIiuWCsqLiMjcmMFwPwxE/TbYNX6cuA12QXwEktMgOR2SUyEl2Cen+/MpacE+A9JCkB6CtGxIC/t9eghSsyFpAX4XMzIEw70w1AdDvf54uB9iw76f8ViwH5n6cVJK8F5Sg/00xxk5kF3o36uIiIiIiCw6hVsiIjLRUC9E66CrFqK1/jhaC9Ea6Krzj0f6Z36NpFQf8iSnwsigD5BigxAbml+fUrN92JWaCUnJgAOXlLC5ifvRAG40xBrqg/jw/O49X+lByBUq9vvsooStECKVULAJsvKvbb9ERERERFaZqw63nHMVwF3ADiAPSDWzt13t64qIyCIwg95WiF6GzhofWEVrJx73t0+6yEG4DHIqoHQPbL3bBzQZER9gZeT4IGf0OCPiR2I5N/X9Y8M+5BrdRgZhZACGemCwxwdSQz1TPx7uB4sHmyUcx4PXTzhOzQwCsWxIy5p0HILULN8mOc2PykpK8cHZ2HHi42Q/kis2lND/xPcxejzoR631tkBPi9/3tkDrWbj0CPS1AzbxzyQzDwo2+y1/kw+8CjZD/kY/ck1ERERERGY073DLOVcMfBh4BZA8ehr/U/vbJrX9Z+DtQI2ZbZrvPUVEZBYGuqDjArSfh/YL/jgxvBoZmNg+LQy5VZBTCZU3+BArJ3gcqYBIuR+BtRCc81MRU9IW5vVWmtiIDw97mv1/j7az0HbO7y88DIf+c2L7cJkPvPI3+LBrbNsA6eGleQ8iIiIiIsuMM7Mrt5p8kXNbgIeAEnyglcjMLHlS+x3AMXzwdZeZPTiv3q4ABw4csMcff3ypuyEiq1k8FoQjtRNDrPbzfutrndg+uwhyq31YlVMVHAfhVW4VZOROPcpKrr2hXv/fcDTwajsH7ef8f9/e5olts4smBl55G4KQsgrCpcH0TRERERGR1cE594SZHZjquTmP3HLOpQLfAkqDU58DPg9sAT461TVmdsI5dwTYDdwNPDjX+4qIrAkjg9DdAF0Nvr5VV73fuusTjhvBYgkXOT/CKn8DbH/R+CifvA0a4bPSpGX7qZ+le57+3GD3xBCz44J/PNWIr6QU///EaJA5GnrlBuFmpHLtjp4TERERkVVnPtMS34YPsgx4h5n9O4BzLusK1z0E7AFumsc9RURWNjPo7wjCqYYp9g0+wOpre/q1qdl+qmCkHDbc4feR8vFAK3cdpGZc+/ck11Z6GMr2+m2y4X7ouOSnOnZeDvbBVNTzD/r/xybU+nL+/6Hc6oRt3fhxTuXCTUUVEREREVlk8wm3fj3Y/3g02JqlY8F+6zzuKSKyfA32+NFU3Q3Bvn7S42A/udYV+Kllo8XaKw/4wCFcNh5eRcogPaJpgzKz1Ewo3u63qYwMBStdBqFX5+Xx7dIv4Mh/jRfiB7/iZKQCirZDyS4/kqxkty90n6yFlkVERERkeZnPT6h78L/+/focrxtdfitvHvcUEbn2RmtbddcHI6umGXE11P30a1OzfTAVLvNF2sOlEC4PzgX7UKmmhsm1kZIWTFfdMPXzsWH///RY6HXJT3lsPuFHfsWHfbvkdB+gleyB0t0++CrZDVn51+ytiIiIiIhMNp9wa/Qn2MZrcC8RkYUVj0Ffuy/O3dMMvS3Bvhl6Wvy+t8Uf9zRNqm0FuOQgqCqDom2w8dnjgdXo+XApZESW5v2JzEdyKuSt89tkI0PQehqajvqt8Sic+T4c/Px4m6wCyFufsG0YP46Uq7i9iIiIiCyq+QROUaAAyJnjdaM/MU9RUObqOefuBj4CJAMfN7MPTnreBc/fA/QBbzazJ2e61jn3t8BLgCHgHPAWM+tcjP6LyAIa6oOOi+MFtxP3nZchPvL0a5JS/RTBUBFkF/vRKKNhVeJUwewifVGXtSUlzY/SKt098XxPMzQegaZjQYH7i1D3JBz/xsTPWFKqr+OVtz4YPbYpYYXHdZCSfi3fjYiIiIisQvMJty7gw60bgU/O4boX4qczHp3HPWfknEsGPgY8D6gFHnPOfdPMjk+6/5Zguwn4F+CmK1z7Q+A9ZjbinPsb4D3Auxa6/yIyT/0dfhTJhC/YF4Li2QnScyB/PZTuhR2/Nh5ShYp9kBUqgoxc1bUSmYtQMWy+y2+JYiPQVRsEzAlb+wWofRwGowmNnV/FcXSFz9GtaJsPwZKSrtnbEREREZGVaz7h1g+BG4DXOOf+wsxarnSBc+75wLPw4db353HPK7kROGtm54P7fRG4F0gMt+4FPmtmBvzSOZfrnCsD1k93rZn9IOH6XwKvWIS+i8iVmPlC2I1HErbDfhTWqOxiX+x603P8lKj8DeP7zDwFVyLXSnLK+JTEyUZXDW0/77e2c+PHx78B/e3jbVOzfU2vsr2+oH3pXijeqZVBRURERORp5hNu/RvwR0AY+Kpz7iVmFp2usXPuLuA/goddzG2012xVADUJj2vxo7Ou1KZiltcCvBX40lX3VESubKgXah6Fiz/3+8YjMNAZPOl8iFV5Axx4a7CK2x4Ilyxlj0VkNpzzxeez8v3qoJONBl/NJ32A3XAYDn0JHvt4cH2yH9VVOhp47fEBWHbhtX0fIiIiIrKszDncMrMa59xfAR8AbgVOO+c+CYzNHXDO3QbsBl4GPBdw+FFbf2RmXQvR8UmmGpJhs2xzxWudc38OjABfmPLmzt0H3AdQXV19pb6KyGQDXXD5l3Dp536rf8rX7HHJftTGrpeNj9wo2Qlp2UvdYxFZDJl5UPEMv/F6fy4eh86LPugaHbV5/kE4/MXx67KLoHiHH9lVvAOKd/kQTAs7iIiIiKwJ81rB0Mz+2jlXDPw+UAT86ehTwf6hhOaj4dFfmdlijNoCP9qqKuFxJVA/yzZpM13rnHsT8GLgrmBK49OY2f3A/QAHDhyYso2IJOjvhEuP+CDr4s/8l1WL+8LTFdfDM38P1t8KVTdBenipeysiSykpabwW166Xjp/vafarNzafgObjfv/k52C4d7xNTrUPu4q2+eL1uet8cfvcakjNvOZvRUREREQWx7zCLQAz+0Pn3IPA+4G9MzQ9hi/K/q353msWHgO2OOc2AHXAa4DXTWrzTeB3g5paNwFRM2twzrVMd22wiuK7gDvMrG8R+y+yug31+pFZFx6CCw9DwyEfZiWn++mFt/8JrLvVH6dlLXVvRWQlCBVD6Dm+zt6oeByil33Q1XQsCL5OwPmfQGxo4vXZxeNB1+iWt86vlBoqUZ0+ERERkRXETTMYaW4v4txefMH49UAO0IMPih4ys8ev+gaz68M9wIeBZOCTZvYB59w7AMzsX51zDvgocDfQB7xltG9TXRucPwukA23BbX5pZu+YqR8HDhywxx+/Jm9ZZPkaGYK6x32Qdf4hqH0M4sN+ZFblDbDxDlj/LD/1SMWhRWSxxePQ0+gXoei8DJ2XEo4vQ2eN/ztqVHbReE2v0fpeBZshKXnp3oOIiIjIGuece8LMpijcukDhloxTuCVrTjwOHReCURLHfQH4y7+E4T7AQfl+2HA7bLgDqm9WvSwRWX7iMehu9MXsm46N1/ZqPjEeeqVk+pp/pXv86K6S3f5xRs7S9l1ERERkjZgp3Jr3tEQRWYP62v0Xv6ZjQa2boM7N8OisXQdF2+G6N/rRWeue6QtEi4gsZ0nJkFPhtw3PGj8/MgStp4OwKwi8jn0dnvj0eJucKr9iY/FOvy/Z5Ud5Jade63chIiIismYp3BKRcbER6G7w03SiNQnTd2r8F7zuhvG2WQX+S9z1bxr/Qle0XTWzRGT1SEmD0t1+47X+nBlEa324Pxr2Nx+Hsz/yq7wCJKdB4TY/sqtwy3hB/LwNkJm7VO9GREREZNVSuCWylsRGoKtuUq2Z0SDrEkTrwGITrwmV+ELLG+4YD7FKdvtiziq4LCJrjXOQW+W3rS8YPz8yCK1ngrArCL0u/BQOf2ni9VkF42HX5C0zT3+vioiIiMzDtOGWcy423XNXycxMoZrIQjODwS7obvIjrKK1Tw+xuiaHVw4i5X5aTdXNsKfKB1k5VZC7DnIqVfBdRGQ2UtITRnklGOqFjou+nlfidukROPxlIKH2aXoO5G+YFHoFj7WCo4iIiMi0ZgqZ9BOUyHIRG/HL23dc8kWPexr9vrsReprGj0f6J13oIFLhA6t1zwxGGyQsex+p9NNuRERkcaRlj496nWx4wP/iof0ctF8YD77qn4Lj35j4y4jULB9yhcsgPXyFLeLbp6T7KZIp6cFxsNeqjyIiIrLKzBRuPcyEXyc+TQ6wP+FxF3Ae6AWygY1AJHjOgINBGxGZylDf+G/3Oy74Lzqj+87LT58umBaGcIn/olPxDAiX+i1U6s/nVPlgS+GViMjylJoBRVv9Nlls2P/dP/rvQPt5aDsHvc3+3GC338YW9JgDl5wQfGX4AC4tC9JC/jh19DgreC7kzyWn+UL5yan+OCklOJcGyYnHweuOhmpjxxkK1kRERGRRTBtumdmd0z3nnNsHfC14+BXgb83ssSna3QD8MfBKfBj2ZjM7cjUdFlnxRgaD1beO+hUHm45By8mJxdrBLy+ftwHK98PuX/fHeev9NMJQCaSHlqL3IiJyLSSnQsEmv80kNgJDPeNh1+g21AOxIf9vTmzQr/w4YT/onx/u9wHZUJ+/pq8dhmv9dMrRLTa4cO8rKWV8BFl62BfYz8hN2Oc9/VyoxI88Tg8vXD9ERERkVZlz7SvnXAHwP0AF8Ptm9k/TtQ0Cr1c7534K/F/gf5xz15tZ+3w7LLJimPmpgk3HghArCLJaTyesqJUOxTtg47OhIFhJK3+D32flL23/RURk+UtO8SHQYq7CGBuB4V4/miw27EOx+Ijfx4b882PHw+Ph2cggjAyM72ND44+H+30IN9AJ/Z3+Fzz9nf5xbGjqfmTk+pArZ3R6fVVQIzI4l5WvumQiIiJr1HwKu78TqAR+OFOwlcjMPuqc+zXgruD698/jviLL02iI1XIyYTsFzSf8D+mjcqp8zZVtLwzqr+zx9VOStb6CiIgsY8kpkJxzbe5l5oOv0dCrvyNYJKUGOmv8vv08XHjIjzRLlB6Bom1QtH18K97up+gr9BIREVnV5vOt+mX4Glr/Pcfrvgo8F/h1FG7JStXT4kdgNZ+AlhM+xGo5CQPR8TaZeVC0w08lLNoOJbuhZKc/LyIiItNzLqj1leWn4U/HzAdfiaFX2zn/b/Lp78FTnxtvmxZOCL22+RUtK2/QNEcREZFVZD7h1rpgP9ephR2TrhdZvkYGfXA1NqXwmN96m8fbZBUEIdYr/NTC0R+cs4v0G2IREZHF5JyfhpiVD2X7nv58b1swkvrE+GjqMz+Ag58Prk+C0j1Q/UxYd4vfh4qu7XsQERGRBTOfcGv0W/uWOV63edL1IstDbxs0HoKGQz7AajwKbWfG62KlZPjQasvzx5dzL96pH4JFRESWq+wCyL4V1t868XxfOzQchEu/gMu/gCc+BY/+i3+uYDNU3wLrnun3eev1yyoREZEVwpnZ3C5w7gngOuA8sNvMBmZxTQZwFNgAHDSzZ8yjryvCgQMH7PHHH1/qbshUzKCr3odYjYf9vuEwdNWOt8mpCqYRBiFWyW7VxRIREVmtRoaCsOsRH3Zd/sV4qYFQKZTt9SO8SnZD6V7/M0FS0pJ2WUREZK1yzj1hZgemem4+39i/ig+3NgBfdc693sw6Z7h5DvAFYCO+Vtd/zeOeInM32AM1j8Kln0P9Uz7M6msLnnRQuMVPRSjd66c0lO7RCoUiIiJrSUoaVN3oN/4A4nE/lfHSI1D7mB/Nfe7H46O5U7N9Hc3SPUHotcc/TsteynchIiKy5s1n5FYWfhTWaO2sNuDTwAPAWaAPyMJPQ3wO8GagED8d8QJ+tFf/1Xd9edLIrSU0GmZd/Jnf6p/0P4y6ZP+DZ9k+KN3n9yW7ID201D0WERGR5W5k0NfvajwSbEf9fnDSYjKhUggVQ7gUQiWT9qUQLlERexERkasw08itOYdbwQtuAX4MVOBHY13xEqAWuMvMzsz5hiuIwq1raLowKykFyq+H9bf5reomBVkiIiKycMz8Co2NR3yx+u4G6G6EnibobvL72ODTr8sq9CPHCzb5Gl8FW/w+fwOkpF/79yEiIrKCLHi4FbxoLvB3wBuAtBmaDgGfBd5lZh0ztFsVFG4tkt5WXydr9LemDYd90XeLK8wSERGR5cUMBjqDoKvR77vrof08tJ2DtrM+ABvlkiC3ejzsKtrqF7Mp2q6SCSIiIoGFrrkFQFBn6+3OuXcDLwZuAMqBENAD1AO/Ar5tZq3zvY+sMbFh6LwMTUd9gDUaZnXXj7fJqfJ1Lna9FKpv9mGWal2IiIjIcuGcn6qYmQfF26duMxAdD7razkLrGb+/9AgM9463C5VA0TYo2uH3xTsUeomIiEwy75Fby41z7m7gI0Ay8HEz++Ck513w/D34umBvNrMnZ7rWOZcPfAlYD1wEXnWl0WcauXUFsWGI1voAa3SL1owfd9X50Vjga2UVbfMF30cLt6rou4iIiKxmZv5npZaTfspjyylf5L7lFAz1jLfLLg5Cr+3BfhsUbvN1v5xbuv6LiIgskkWZlricOOeSgdPA8/C1vR4DXmtmxxPa3AO8Ex9u3QR8xMxumula59yHgHYz+2AwQi3PzN41U1/WfLg1EA3CqxofWkVrE7YaX5NiNLwCPww/UuGH4udU+X3eOl/wvWgHpGYs3XsRERERWS4SQ6+Wk9Ac7FtPw2DXeLuM3CDwCqY2FgajvSLlCr1ERGRFW5RpicvMjcBZMzsP4Jz7InAvcDyhzb3AZ82neb90zuU658rwo7Kmu/Ze4M7g+s8ADwIzhlurxS8++xekNx+aVdtw0hDrUztI7amf+MMVQFIq5FT44GrD7T68StwiFZCcugjvQERERGQVcQ5yq/y25Xnj5838Lw9bTvmtNdif/A48+dnxdlkFwSj4vcEK0nt9Yfuk5Gv/XkREZF4GhkY4efwgXUe/R1bTkyTZyKyuy3zW77Djphcscu+W1moJtyqAmoTHtfjRWVdqU3GFa0vMrAHAzBqcc8VT3dw5dx9wH0B1dfU838LykhS9TEHf+Vm17Y6l8KAVEi55Abt27iJcvN6HWTmVfsh8UtLidlZERERkrXLOj8qKlMOmZ098rrctGOV1HBoO+cV5Hv1XiA3551Oz/Gj50RIQJbv9yo1ZBRrlJSKyDLR0D3LwXA3tRx8gXPsQu/sfY79rBqDJFTLgMmf1Op29q35tv7mHW8652SUe0zMz23SVrzHZVP/6Tp5vOV2b2Vw7IzO7H7gf/LTEuVy7XN30zs/Mum1dZz+ff+AM//VELWmNSbzpmet5xx0byc2aaRFNEREREVlU2QWQfSusv3X8XGzYj+xqPDy+eM+Rr8DjnxhvkxaC3HWQt/7pW261ykaIiCyCkVic0009HLzcQd3pxwnVPMTegce5I+kkaS5Gv8ugNu8AJzb+NmXXv4iSym2zfu11i9jv5WI+I7fWM30olGg05JncbjHCn1qgKuFxJX61xtm0SZvh2ibnXFkwaqsMaF7QXq8SFbmZfPDle/mtOzbx4R+d5t8ePscXfnmJtz9rI2+9bT3hDE07FBEREVkWklOhdLff9r/OnzODjot+lFfHJX/ccRHaz8O5H8NI/8TXCJX4wvXZxZBdBKEiv88unnicXajyEyIiUzAzLrb1cbi2k0M1UQ7XdBBrOMSL7WHuSX6UMtcOQHtkM+0b3kb+vnvI3PBMtqSkL3HPl685F5R3zl3kygFVEpAPZAePDWgAhgHMbMOcbnrlPqXgi8LfBdThi8K/zsyOJbR5EfC7jBeU/79mduNM1zrn/hZoSygon29mfzpTX9Z8QXngZGMX//jD03z/WBN5Wan89p2b+I1b1pORqpoOIiIiIiuKGfQ0+7Cr89L4vrfVn+9thd5mGBmY+vrMvIlhV3bRpEAs2CLlkDq76TUiIivJcCxOXUc/p5q6OVzbyeHaKIdro0T7hymig1emPcKr037OupGLxJJSGVh/F1m77sFtfq6vXy1jlmy1ROfcbuD3gbcBDwMvN7O2RbrXPcCHgWTgk2b2AefcOwDM7F+dcw74KHA30Ae8xcwen+7a4HwB8GWgGrgMvNLM2mfqh8KtcYdrO/m7H5zm4dMtFIfT+d3nbOZVB6oUcomIiIisJmYw1BOEXS1+GzsOwq/e1vHzA51Tv05Woa/ZmlM5Xr81p9IX0c+p8s+rlquILEN9QyNcbu/jYmsfl9t7udTW5x+39VLfOUAs7nOX5CTHnuI0XhU5wp39P6Ks9RGcxaHiAOx7Dex+OWTlL/G7Wb6WLNxK6MDb8TWpfgHcbmaxRb/pElG49XSPnm/j739wml9dbKc4nM59t2/kdTdVk5W2WtYzEBEREZFZGxmCvrYg/GqGnhboqoNoDURr/dZZA8O9E69LTvdB11gNsA0Jx+sgPXzt34uIrBnxuFHX2c+5lh7ONvdwrqWXcy09XGjtpaV7cELb3KxU1uVnUV2Q7ff5mey1U2yu/yYpJ74Bg1GIVMDeV8O+10LR1iV6VyvLkodbQSceAO4E/peZ/ds1uekSULg1NTPjF+fa+OhPzvLIuTbyslJ5220beOMt68nJVC0GEREREUlg5kd4jYZd0VrovDw+NbL9ov9ymCircGLh+9ERX6MjwBR+icgs9A/FuNDqgyu/9XK2uYcLrT0MDMfH2uVmpbK5KMSGwmzWF2ZTnZ/FuoIs1uVnk5OVCvEY1DwKJ/7Hb9Eav0rtjl+D/a+F9bdrNOocLZdw63eB/ws8Yma3XZObLgGFW1f2xKUOPvaTs/z4ZDPh9BTe9Mz1vPW2DeRna3VFEREREZml/o7x4veJW/sFPxIsPjKxfUbOxLBrdPpj/gY/CkxTgUTWDDOjsWuA88Hoq8R9Xef4IhrOQWVeJpuKQmwuCrGpOOSPi0NTf3+NDcOFh32YdfLbfnRqcjpseg7svBd2vATSQ9fwna4uyyXcegW+flW7mRVek5suAYVbs3esPso//+Qc3znaQEZKMq+7qZr7bt9ISUTLS4uIiIjIVYjHoKcpGPU1abrj6LnJtb8ycseDrsn7cJlGWIisUL2DIxyr7+JwbSdH66KcDUKsvqHxaknZaclsLAqxsSibjYV+vyl4fMWa0cP9fmXZE/8Dp74DA1FIzYatz/dh1pbna+ToAlku4dbv4Yu295tZ9hWar1gKt+bubHM3//zgOb5xsJ5k53jFgUped2M1u8oj+HUAREREREQW2GC3n+rYfgE6Lkzcd16GxDLBKZm+Jk7RDijePr7PqVboJbKM9A2NcLy+iyN1UY7URjlcF+VcSw+jsUdpJIOtpWE2FWWzsSjEpkK/L4mkz/67Z28r1D0J9U/6/cWf+RqBGbmw7R4faG16tlaAXQRLHm4559KAx4A9wBkz27boN10iCrfm73JbH//68Dm+8ngtQ7E4m4tDvHR/Offur6AqP2upuyciIiIia0VsxI/uGg272s5BywloPgnd9ePtUrPHQ6+ibVC4FULFkFUA2YWQFvLzmkRkwQ2OxDjZ0M3h2k4O1fow60xzN8HChBSH09lbmcOeilz2VuawuyKHonD63G4yEIX6g+NBVv1T/u8GAJz/3K97pg+01j8LklVPejEtWbjlnEsGngX8H+BWwIB/MrM/WLSbLjGFW1evs2+I7xxp5OtP1fGri+0A3LA+j5deV8GL9pSRm6XaXCIiIiKyRPo7oeXUeNjVEmzdDU9vm5weBF0FvuD9aOiVVQDJaf6LcFIqJKcE+9SEc8GWmQ+hEsgughT9HCxrUzxunG/t4VBNlEO1nRyq6eREQzdDMV/gvSA7zQdZlbnsrchhT2XO3Mvd9LVD4xFoPAwNh32Q1XZm/Pm89VB+PVRcD+XXQdk+TTe8xhY03HLOnZ9l0zSgEBiNLh3QDuwxsyn+5l8dFG4trJr2Pr55qJ6vPVXH2eYeUpMdd24r5mXXVfCc7cVXnv8sIiIiInIt9HdA23noa/XTlvpaoa8NetsmnWuHwa753WM06AoVB1vJ+D5S7ovkRyogZY6jU0SWETPjcnsfR+v89MJDNb5WVvegXyQiOy2ZPZU57KvKZV9lLvuqcinPyZj9tEIzX3svMchqPJwwIgv/OSrbF4RZ1/m9Fp1YcgsdbsXxI7DmOr72LPAaM3tyjtetKAq3FoeZcay+i68/Vcc3DtXT0j1IOD2FO7YV8Zztxdy5rVirLYqIiIjIyhAbhtiQ38dHJh0PQ3x4vE1fuy+O39Mc7IPj3mboboKR/qe/fqhk4oqQk1eJzCrQdElZFmJx40JrL8fq/bTCo/VRjtV30T3gg6zUZMeOsgj7Kv3Uwv1VuWwsCpGcNMv/f+NxaD8HDYfGt8Yj0N8eNHBQuAVK90DpXijb6/fZq3YNvBVtocOti/hw60oGgU7gOPA94OtmNjSnm61ACrcWXyxuPHKulW8erOcnp5pp7RnCOdhflctzthXz7O3FKkYvIiIiIqufGQz1+JCru37iapCjK0RGa58egKVk+JEpORU++Bo7roRIpT/WdCtZYPG4cbGtl4M1nRyujXK0Lsrxhq6xVQvTU5LYURZhd0WE3eW+RtaWkhDpKbOcrRMb9lOEGw6PB1lNR/1nBPxU4OIdfkRW6V6/L9kFaat2vbtVZ8kLyq8lCreurXjcOFIX5ccnm3nwVDOHaqMAlETSefY2P6Lrti2FhNJTlrinIiIiIiJLwMyP/orWBFud33fVBce10NMIFp94XVrYT8PKLvT1wrIL/eOx46COWEYE4jG/umQ85l/H4pPOxXw/UtJ9HbKUxC3Dhw6je60+uWq09QxysKZzbDtU00lXMCIrKy2ZXeURdgUh1u6KCJuLQqQkz+G/f2cN1DwKNb+C2l9B03GIDfrnUrP9aKyyfX40Vtk+KNymunUrnMKta0jh1tJq6R7kwVPN/ORUMz893Ur34AhpyUncua2I19xYxe1biub2F6aIiIiIyGoXG4buRh90dQWBV3ejrxk2Vi+s3R+PDCxuX5LTgpArZfx4tLj+hCL8QRCWlAIuGZKSg33SpMfJ4+HZaJj2tH1wnJoJqVl+n5YdPB7dZ2oq5wwGR2IcreviqcsdY2FWbYcfMZjkYFtphP1VueyvymF/VR6bi+cwtRD8/6ONh32QNRpoddX551KzoOIZUL4fyvb7ICt/o/9vL6uKwq1rSOHW8jEci/P4xQ5+dKKJbxyso7VniNJIBq94RiWvOlBFdUHWUndRRERERGTlMIOh3oTQq80Xx08Mksb2buI553wNsZHBhG3Aj7SZfG6sDtkQxBKPR+uRjR7HfNvEUWOTR4zF48F9B8Zff1ZVdqaQGHyl50BGjh+5lh6ZYh88n5nnR7xl5vupnqskIGvpHuTJyx08cclvR2qjYysXVuRmsq8qJwiz8thdESErbZYzaWIjPrTqvAQdl6DtLNQ+BnVPjk+vzamGqhuh6ia/L9ntVxyVVW8xVks04LfM7EdzuO524NOAmdmmOd10BVG4tTwNx+I8cKKZLz12mYdOtxA3uHVzAa86UMULdpVq1UURERERkbXAzAdjiWHXyKAPToYHYLhvfBvqg+F+GO71+6He4HwvDHTBQNSHewNdMBiFwe6nT+9MlJTqw67EwCsrz++zixK2Qr8KZlbhsphGF4sbp5u6eeJSB09e6uCJyx1causDIC05iT2VOTxjXR7XV+dxfXUuxZGMp7+Imf+zG+jyf059bT7A6rzsQ6zOYIvW+VByVFKqH4k1GmRV3ehXBpU1abFWS3yZmX1zDte9APguPtxatUmCwq3lryHaz1cer+VLj9dQ29FPTmYqL7uuglffUMWOsshSd09ERERERFaieNwXLx8NvAaiflW+/g4/rXPCccf4cV/beK2oyTJyxkOvrAI/emzK6ZUJx8npwXRNBy7p6VtScnAcjK4bm+qZAkmpDJPMmdYBDtb18kRtD0/UdjMwMEi2G6Aya4TripPYWZDE1lyoyBohdaTXB1aDwXsf7A5CrK6E8K97YmiVKFwGudWQuw7y1k08jlT4/omgcOuaUri1csTjxi/Ot/HFx2r4/tFGhmJxblifx+88ezN3bC3SaosiIiIiIrL4Rle97G2Bnha/723xtc56mxOOW/0IswkjzgZmHi12rSSn+2mX6aHxaZnpEX9udKpm4nFmrg+wcqogdYqRXiJTmCncupYTUzOD/TSRtMi1lZTkuHVzIbduLqSjd4j/fqqOT/z0PG/+1GPsqcjhd569mefvLCFpLoUORURERERE5sK5IBgK+0LocxUbGQ+7YqN1y+LjK1eObb4uWd/gMKcbo5yo7+RUfScXmjux2DApxFifm8a24gy2FmayqSCdSFowjTMpOQioIkGAFfQ3PQJpoWUxfVLWtmsZbt0c7Fuu4T1FZiUvO4233baBN968jq89Vcu/PHiOd3z+CbaWhPidZ2/mRXvKtMqiiIiIiIgsP8kpkBzyodMUmroGeOxiO49f7ODxS+0cr+8ibpDksthZXspNNxVw04Z8blifT162QipZmWacluic2wvsn3T60/hpif8EPHml1weygeuBNwCpwNfN7OXz6+6UfcwHvgSsBy4CrzKzjina3Q18BEgGPm5mH5zpeufc84APAmnAEPAnZvbjK/VH0xJXh5FYnG8faeBjPznL6aYe1hVk8b/u3MTLrqskLUUhl4iIiIiILD8jsTjnWnp54lIHj19s57FL7dS0+1UGM1OTua46lwPr87lhfR7XVecRStcqg7JyzLvmlnPuL4H3Tj4d7Oe6fqoLrrnbzH44x2unf1HnPgS0m9kHnXPvBvLM7F2T2iQDp4HnAbXAY8Brzez4dNc7564Dmsys3jm3G/i+mVVcqT8Kt1aXeNz4wfEmPvaTsxypi1Kek8Fv3bGJV99QpRUWRURERERkyXQPDHOysZvj9V0cr+/iRGMXJxu7GRrxNbgKQ+ncsD5vLMzaURYhVbNRZAW72nDrLxeoH83AX5jZvy/Q6wHgnDsF3GlmDc65MuBBM9s2qc0twPvM7AXB4/cAmNlfz/J6B7QC5WY2Y80whVurk5nx0OkWPvaTszx2sYOC7DTe9Mz1vOHmdeRr6K6IiIiIiCwSM6Oxa2AsxDre4LdLbX1jbfKyUtlVnsOOsjA7yyNcV5XHuoIsLZIlq8rVFJT/On6qXqJP4UdgfZQrT0uMAz3ABeCI2XRrf16VEjNrAAgCquIp2lQANQmPa4Gb5nD9y4GnrhRsyerlnOPObcXcua2YX55v498eOsc//PA0//zgWV51oIq33baBdQXZS91NERERERFZwWJx40JrL8fqoxyv7+JYEGa19w6NtdlQmM2u8givfEYlO8sj7CzLoSSSriBL1rQZwy0zOwQcSjznnPtUcPiAmX1zsTo26Z4/AkqneOrPZ/sSU5yb1bRK59wu4G+A58/Q5j7gPoDq6upZdklWqps3FnDzxgJON3Xz8Z+e54u/quFzv7zE3btKue/2jVxXnbfUXRQRERERkWVuJBbnVFM3h2ujHKuPcqy+i5MN3fQP+zEhaclJbC0N8bwdJeyqiLCzLML2sojqZIlMYT6fircE+yuN2lowZvbc6Z5zzjU558oSphU2T9GsFqhKeFwJ1AfH017vnKsEvgb8hpmdm6F/9wP3g5+WONv3JSvb1pIwH3rFPv74+dv49CMX+fwvL/Hdo43cuD6f37x9I3dtLyYpSb89ERERERFZ68yMmvZ+DtZ2cqjGb0frowwM+/pY4fQUdpRHeM2NVewqz2FnWYTNxSEtZiUySzPW3FoJnHN/C7QlFITPN7M/ndQmBV9Q/i6gDl9Q/nVmdmy6651zucBDwF+Z2Vdn2x/V3Fq7egZH+PJjNXziZxeo6+xnY1E2b79tI79+fYWKz4uIiIiIrCFtPYMcrouOBVmHaqNjUwvTU5LYXZHDvspc9lX5fXV+ln4xLnIF8y4ovxI45wqALwPVwGXglWbW7pwrBz5uZvcE7e4BPgwkA580sw9c4fr/DbwHOJNwu+eb2VQjw8Yo3JKRWJzvHG3k/ofPcbSui/zsNN5w8zreePM6isLpS909ERERERFZQK09gxypi3K0Nur3dVHqowMAOAdbi8M+xKrKZV9lLttKw1q1UGQeVnW4tdwo3JJRZsYvz7fziZ+d50cnmklLSeJl+yt427M2sLUkvNTdExERERGRORoNso4kBFkNQZAFsLEwm90VOeypyPH7yhzVyBJZIPNaLdE5N7qyoZlZyhTn52vC64msVs45btlUwC2bCjjX0sMnf3aBrz5Zy5cer+GOrUW8/VkbuG1zoVY1ERERERFZZsyMpq5BjtZFOVrvQ6yjdV00dk0Msm5Yn8+eIMTaVR4hnJG6hL0WWbumHbnlnIsHh2ZmyVOcn68Jr7faaOSWzKS9d4j/ePQSn/nFJVq6B9lWEuZtz9rAvfvLSU9ZtR8LEREREZFlKx436jr7OVbfxbH60RFZXbT2DAJ+auGmohC7yyPsDkZkKcgSufbmNS3ROfcgYABm9uypzs9X4uutNgq3ZDYGR2J882A9n/jZBU42dlMUTue3bt/I629aR2aaQi4RERERkYUWixs17X2cae7hTHM3Z5t6ON3czdnmnrFVC5OTHFuKQ+wqz2FPhQ+zdpRFyNbUQpElp5pb15DCLZkLM+PnZ9v45wfP8si5NgpD6bzjDoVcIiIiIiJXo6V7kKP1UY7Xd3G6qZszTT2ca+lhcGR8IlJZTgabi0NsKQ6ztSTE1tIwO8siWulcZJlSuHUNKdyS+Xr0fBsfeeCMQi4RERERkVkyM2o7+jlWH+VYfRdH6/y+uXtwrE1FbiZbSkJsKQ6xpSTMluIQm4tDmlYossIo3LqGFG7J1frVhXY+8sBpfn5WIZeIiIiIyKjRIOtQbSeHajo5WudrZHUNjAB+SuHmohC7KiLsKvd1sXaWR4goxBJZFRRuXUMKt2ShTAy50njHHZsUcomIiIjImtHeOzQWZB2q6eRQbZT23iEA0lKS2FEWYVe533aX57CtNKwphSKrmMKta0jhliy0ySHXO5+zhdfeWE1aStJSd01EREREZEF09A5xoqGL4w1dHKqNcqimk8vtfYBfrXBLcYh9lbnsq8plf1Uu20rDpCbr52GRtWS+qyW+d7E6ZGZ/tVivvdQUbslieexiO3/3/VM8eqGddQVZ/PHzt/GiPWUkJbml7pqIiIiIyKyMxOJcaO3leEMXJxq6OdnYxYmGLpq6xmtkledksK/KB1n7KnPZU5lDSKsViqx58w234sCiDOsys1U7VlThliwmM+PB0y38zXdPcrKxmz0VObz7hdu5dXPhUndNRERERGSMmVHX2c+Z5h7ONvVwqskHWaebehgKVixMTXZsKgqxoyzCjrIwO8oibC+NUBROX+Lei8hydDXh1mIwhVsiVycWN75xsI6//8Fp6jr7edaWQt5193Z2V+QsdddEREREZA2Jx0dDrG5ON/VwpqmHs83dnG3uoXcoNtauMJQWhFfhIMyKsKkopFIbIjJr8w237lisDpnZQ4v12ktN4ZZcSwPDMT7/y0t89Cdn6ewb5qX7y/mj52+jKj9rqbsmIiIiIquMmXGprY9DtZ0cDIq8n2jopn94PMQqiaSzpTjM5uIQW0pCY8f52WlL2HMRWQ1UUP4aUrglSyHaP8y/PXSOT/78ArG48Yab1/HO52zRDxEiIiIiMm/N3QMcromOhVmHa6NE+4cByExNZndFhN0VOWwrCbOlJMTmojA5WalL3GsRWa0Ubl1DCrdkKTVGB/jIA6f50mM1ZKel8I47N/HWWzeQmbZqZwKLiIiIyFUarY91vN6vVni8voujdVHqowMAJCc5tpWE2VeVM7Zi4ZbiEClarVBEriGFW9eQwi1ZDs42d/M33zvFD483URJJ5/953lZe8YwqkrWyooiIiMiaNjQS50xz91iQdSIIs7oGRgBwDjYUZrOrPId9lTnsr8plV3mOflkqIkvumoRbzrkwUAqEgB6g0cy6F+TFVxCFW7KcPHaxnf/vOyd46nInW0tCvOvu7TxnezHOKeQSERERWe3ae4c4kRBgHW/o4lxLD8Mx/x0wMzWZ7WVhdpZF2FkeYWdZhG2lYbLSUpa45yIiT7do4ZZzrhL4beBlwFYg8RuzAaeBrwL/amZ1877RCqJwS5YbM+N7Rxv50PdPcaG1lxs35PNn9+xgf1XuUndNRERERBZALG5caO3heEP3WJh1oqGLpq7BsTYlkfSxVQp3BUHWuoJsjewXkRVjUcIt59z/Av4GGF2Wbaq/FUdfvA/4UzP7l3ndbAVRuCXL1XAszhd/dZmPPHCG1p4hXrSnjD95wTbWF2YvdddEREREZJb6hkY42einFR4LRmOdauxiYDgOQGqyY1NRaGw01migpYWGRGSlW/Bwyzn3XuAvRx8CMeAEcBboBbKBzcAOYHRytgHvM7P/M+cbztyXfOBLwHrgIvAqM+uYot3dwEeC/nzczD44m+udc9XA8aDvf3el/ijckuWuZ3CE+x8+z78/fJ7hWJzX3ljN7z5nMyWRjKXumoiIiIgkaO0ZnBBiHa+PcqG1l3jwFS4nM3XClMIdZRE2F4dIS1GhdxFZfRY03HLO3Qz8DEgCRoC/Bz5sZk1TtC0Bfh/4IyAViAO3mtmjc7rpzP35ENBuZh90zr0byDOzd01qk4yfIvk8oBZ4DHitmR2/0vXOua8G/X5U4ZasJs1dA3zkgTN86bEakpMcv3HLOt5xxyYKQulL3TURERGRNWVgOMbZ5h5ONHRxqrGbk43dnGzsorVnaKxNZV7mhCBrV0UO5TkZqqUqImvGQodb/wG8Bh/4vNzMvjGLa14CfD14+EUze/2cbjrza58C7jSzBudcGfCgmW2b1OYW/MirFwSP3wNgZn890/XOuZcCt+JHo/Uo3JLV6HJbHx954Axfe6qWjNRk3nrrBn7zWRvJyUpd6q6JiIiIrDrNXQMcqo1ysqGLk03dnGzomjAaKyM1ia0lYbaXhtlW6oOsnWUR/WwmImveQodbNUA58BUze/Ucrvsi8Cqgzsyq5nTTmV+308xyEx53mFnepDavAO42s7cHj98I3GRmvzvd9c65bOBH+NFef8wM4ZZz7j7gPoDq6upnXLp0aaHensg1c7a5hw//6DTfOtxAOCOF+561kbfctoFQulbLEREREZmPvqERjtRGOVjTyaHaTg5e7qQ+OjD2/LqCLLaVhNleFmFHaZhtpWEVeRcRmcZM4dZ8vrUWBfvvzfG67+PDrcK53tA59yOgdIqn/ny2LzHFuSuleu8H/tHMeq401NfM7gfuBz9ya5Z9EllWNheH+Ojrrud/3dnFP/zwNH//w9N86pGL/PYdm3jjLevISE2+8ouIiIiIrFHDsTjnWno4eNkHWU9d7uR0U/fYiKzq/CwOrM9nX1Uu+6ty2F4aIVu/RBQRWRDz+du0BT9yq2+O1422b53rDc3sudM955xrcs6VJUwrbJ6iWS2QOFqsEqgPjqe7/ibgFUFNrlwg7pwbMLOPzrX/IivJzvIIH3/TAQ7WdPL3PzjFB75zgn//6Xl+645NvPqGKo3kEhERkTVtOBbnUlsvp5t6ON3UzZlgf6G1l5EgycrJTGVfVS7P31XKdVW57K3MUV1TEZFFNJ9vqU/iw609+FUGZ2tPsH9iHvecyTeBNwEfDPZT1QB7DNjinNsA1OFrhr1upuvN7FmjFzvn3oeflqhgS9aM/VW5fO5tN/HL8238ww9P83++dZwP//A0r7mxijffuoGK3Myl7qKIiIjIojEz6jr7OVbfxcmGbk43d3MmCLGGYz7Ecs6PyNpSHOa5O0vYVhJmX1Uu6wuyVOhdROQamk/NrXuAb+FHcO00s7ZZXFMIHMNPSXyRmc11SuNMr10AfBmoBi4DrzSzdudcOfBxM7snod8fBpKBT5rZB2a6ftI93ocKyssa99TlDj7xswt892gjAC/cXcrbn7WR/VW5S9sxERERkas0EotzrqWXY/VRjtd3cay+i+MNXUT7hwEfYlXlZbG1JMTm4jBbS0JsLQmzqShEZppKN4iIXAsLWlA+eMGPAb8NHAFebWYnZ2i7DT/Cay/wUTP7vTnfcAVRuCWrXV1nP5955CL/+ehlugdHeMa6PN5+2waev6tUxU9FRERk2esfinGisYtjddGxEOtkYzdDI3EA0lOS2B6sULirPMLO8gjbS8Nkpak0g4jIUlro1RJvDw7/AHgpMIwvFv8AcBZfWysL2Aw8B7gbP/3x68BHZnptM3t4Tp1ZhhRuyVrRMzjCfz1ewyd/foGa9n4q8zJ5y60beNWBSsIZWqpaREREll60b5hj9T7EOlYf5Wh9F+dbesaKvOdmpfoAqyzCrvIcdpVH2FCYTUpy0tJ2XEREnmahw604E1cadMy88uCVnh9lZrbifx2icEvWmljc+OHxJj7xs/M8drGDcHoKr7qhijc/cz1V+VlL3T0RERFZI1p7BjlSF+VYXZSjdV0crY9S29E/9nxZTga7ysdDrF0VOZTnZKg2lojICrEY4dZiMDNb8RPWFW7JWnaoppNP/vwC3z7cQNyM5+0s4a23buDGDfn6wVFEREQWTHPXAEfroxyp7eJIXZSjdVEauwbGnl9fkMWuihx2jwZZ5RGtVigissItdLj1lwvSqymY2fsX67WvFYVbItAYHeBzv7zIFx69TGffMLvKI7z11g28eF8Z6SkrPsMWERGRa2Q4FudCay+nm7o53djNsXofZjV3DwK+0PvGwmx2V+SwpyKH3RU57CyPEFGJBBGRVWfBC8rL9BRuiYzrH4rx9YN1fPJnFzjT3ENhKJ033ryO199cTaF+eyoiIiKBWNy41NbL6aYeH2QF24XWXoZj/vtKkoNNRSH2VOSwKwizdpZHCKWv+MomIiIyCwq3riGFWyJPZ2b87Gwrn/jZBR481UJaShL37ivnN2/fyNaS8FJ3T0RERK6haP8wJxq6OB6sVHiioYuzzT0MjoxXP6nKz2RbSZgtJeFgH2JTUYiMVI0AFxFZqxRuXUMKt0Rmdra5h08/coGvPFHLwHCcZ28r4r7bN3HzRtXlEhERWU3MjProgA+xgtUKjzd0TSjyXhhKY0dZhO2lYbaW+G1zcYhsjcYSEZFJFG5dQwq3RGanvXeIz//yEp955CJtvUPsrczhvts3cveuUi2/LSIisgJ19g3xVE0nBy938lRNJ4dqOon2DwO+NtaGgmx2BMXdd5ZF2FkeoTicscS9FhGRlWJRwy3nXBKwCcgDZvWvk5k9fFU3XcYUbonMzcBwjK8+WcvHf3qBC629VOVn8rZbN/CqG6rIStNvbUVERJaj4VicU43dPHW5g6cud3KwppPzrb2Ar421tSTM/qpcdlXksDMYmaXRWCIicjUWJdxyzt0F/CFwF5A2h0vNzFbtv2wKt0TmJx43fniiifsfPs8TlzrIzUrljTev4zduWU9RWMXnRUREloKZ0RAd4GxzD+da/HayoZsjddGxGlmFoXSuq87luupc9lflsrcyV0XeRURkwS14uOWc+xDwR6MP53i5mdmqrQSpcEvk6j1xqZ37Hz7PD443kZqUxIv2lvGGm9dxfXWu6nKJiIgsguFYnIutvWMhlt/3cq6lh76h2Fi7SEYKW0vC7KvyQdZ11blU5Gbq32cREVl0CxpuOedeBXwx4dQZ4GdAEzA4m9cws/fP6aYriMItkYVzvqWHTz9ykf9+so6ewRF2lkV44y3ruHd/uaYsioiIzFPXwDAnG7o5HhR4P97QxenGHoZi46sVVuRmsqk4xKaibDYVhdhc7FcrLAylKcgSEZElsdDh1kPAs4Bh4G1m9vmr7+LqoXBLZOH1Do7wjYP1fPYXFznZ2E04PYWXP6OSN9xczebi8FJ3T0REZNlq6hrgSG0QYtX7IOtye9/Y8wXZaews98Xdd5RG2FwcYmNRtn6JJCIiy85Ch1udQBj4FzP73avv3uqicEtk8ZgZT17u4HO/uMR3jjQyFItzy8YC3njLOp63s4RUrbIoIiJr2MBwjCN10WC1Ql/ovSE6AExcrXB0pcJdZRGKwukaiSUiIivCQodbUSAEvM7MvrQA/VtVFG6JXBttPYN8+fFavvDoJWo7+ikOp3Pv/nLu3V/BrvKIflAXEZFVzcy40NrLwZpOngrCrJMN3YzE/c/2VfmZXFeVx/6qXPZV5bC9NKLVCkVEZEVb6HDrELAbeLOZfW4B+reqKNwSubZiceOh0838x6M1PHS6meGYsbk4xL37fNBVXZC11F0UERGZt3jcqOvs52xzD2eauznT1MOZZl/wvWdwBIBQegp7K3P8ioVVeeyvzqUwpJWGRURkdVnocOv/AH8O3G9m71iA/q0qCrdElk5n3xDfOdLI1w/W8asL7QBcX53LvfsreNHeMv2gLyIiy5aZUR8d4Hh9F6ebusfCrHPNvfQPj69WWBROZ0txiC3FIXaURbiuOo/NxSGSkzRiWUREVreFDrdKgCP4qYk3mNmxq+/i/Dnn8oEvAeuBi8CrzKxjinZ3Ax8BkoGPm9kHr3S9c24v8G9ABIjj3+/ATP1RuCWyPNR39vPNQ/V8/ak6TjZ2k5zkeNaWQu7dX85zd5QQzkhd6i6KiMgaNRyLc66lxxd4r+/iWFDoPdo/PNamPCeDzSXhsSBrc7DlZqUtYc9FRESWzoKGW8EL3gJ8H+gDfsfMvnp1XZw/59yHgHYz+6Bz7t1Anpm9a1KbZOA08DygFngMeK2ZHZ/ueudcCvAk8EYzO+ScKwA6zSzGDBRuiSw/pxq7+cbBOr5xsJ66zn7SUpK4fUshL9xdxnN3lJCTpaBLREQWR/9QjBONXRyti44FWaeauhkaiQOQnpLE9rLxIu87yyJsLQnplzAiIiKTLHi4FbzoOuDrwF6gCXgCaMOPcJqJmdnb5nXTqftxCrjTzBqcc2XAg2a2bVKbW4D3mdkLgsfvCTry19Nd75y7B180/w1z6Y/CLZHlKx73qy1+92gj3z3SQH10gJQkxzM3F3LP7lKet7OEAk1dFBGReeofinG8IcqR2ihH6nygdbalh1hQ5D0vK5Vd5TljIdau8ggbCrNJ0Wq/IiIiV7QYI7dygA8DrwdGl12Z9QuZWfKcbzp9XzrNLDfhcYeZ5U1q8wrgbjN7e/D4jcBNZva7013vnPsD4BlAMVAEfNHMPnSl/ijcElkZzIxDtVG+e7SB7x5p5HJ7H0kObt5YwAt3l/KCXaUURzKWupsiIrJMRfuGOdnopxMeqYv6IKu5hyDHojCUxu6KHPZU5Izty3IytJqviIjIPM0Ubs15PWDnXAj4MbB/8lOzfIk5p2nOuR8BpVM89eezfYl59CMFuA24AT/98oHgD/KBKfp3H3AfQHV19Sy7JCJLyTnH/qpc9lfl8u67t3O8oYvvHmnkO0cb+ItvHOO93zzGvspc7txWxB1bi9hbmativSIia9BILM7Ftl5ONHRzsrHL7xu6qI+Ol2EtDKWzpyLC3btKfZBVmUNpREGWiIjItTLncAv4PeC64Lge+Cjwc/zUxMEF6tcEZvbc6Z5zzjU558oSphU2T9GsFqhKeFyJ7zvAdNfXAg+ZWWtwn+8A1wNPC7fM7H7gfvAjt+b27kRkqTnn2FWew67yHP7o+Vs509zDd4808pNTzXzkgTN8+EdnyMtK5VlbirhzWxG3by3SyosiIqvQwHDMj8SqjXKsPsqJhm5ON3UzGNTHSklybCoKccOGfLaXRtheFmZHaYSSSLqCLBERkSU0n9USDwO78SsL3jga/iwV59zfAm0JBeHzzexPJ7VJwReUvwuowxeUf52ZHZvueudcHj7Iug0YAr4H/KOZfXum/mhaosjq0t47xE/PtPDQqRYeOt1CW+8QAHsqcsZGde2vylW9FBGRFWZgOMaJYEqhr5EV5UzzeH2sguw0dpRF2FEWHguyNheHSE9ZsOoaIiIiMgcLWnPLOdcDZAJ/amZ/vwD9uyrBKoZfBqqBy8ArzazdOVcOfNzM7gna3YOvE5YMfNLMPjDT9cFzbwDeg5/C+J3JodlUFG6JrF7xuHGsvosHTzXz0OkWnrzcQdwgMzWZPZU5XFeVy3XVueyvyqM0R/W6RESWAzOjPjrA2eYezjb3cKqxi8O1Tw+ydlfksLcyZ2yvaYUiIiLLy0KHW01AIfBqM/vKAvRvVVG4JbJ2RPuG+enZFh6/2MFTNZ0cr48yHPN/p5ZGMtg/Fnblsqcyh6y0+cwEFxGR2RiOxbnU1sfZ5h7OtfSMhVnnWnroG4qNtcsfDbIqxoMsFXoXERFZ/ha0oDxwHLidqQu8i4isGTlZqbx4bzkv3lsOwOBIjOP1XRys6eSpy50crOnke8caAUhOcqwvyKI8N5OynAzKcjIpz524z05X+CUiMhsjsThnmns4VNPJodpODtZEOdvcPfYLBoDynAw2FYd41YEqNheHxraC7DQFWSIiIqvMfEZuvRX4OPBTM7tjUXq1gmnklogkausZ9F+8LndyuqmHhmg/9dEBWnsGmfzXbyQjhfLcTIrC6WSkJvstJSk4Tho7lx6cy05PpjCUTnE4g+JwOrlZqfrCJiKrjplR29HPodpOH2bV+PpY/cN+NFZOZir7qnLZVR5hSxBgbSwKEdIvDERERFaVhZ6W6IAfAs8G/szM/ubqu7h6KNwSkdkYGonT1DVAQ3TAB16d4/vWnkEGhmMMjsQZGI4FW5yBkdjTArFEaclJFIXTKQqnUxxOpzgyHnwVR9IpCmVQFE6nIJRGqgrgi8gyZGY0dg1wpDbK0TofYh2ujY4t5pGWksTu8gj7qvyU732VuawryFKwLyIisgYs6LREMzPn3L3Ap4D/zzl3O/Ax4FEza7u6roqIrA1pKUlU5WdRlZ8162vMjKFYnIHhOIPDMXoGR2jpHqR5bBugpcsfX2zr5VcX2+nsG57ytfKyUseCsKKQ3xeG0slKTyE9JSnYkklPTSI9OcnvU5LHzoczUsjJTCUpSV8oRWR+Rgu9JwZZx+qjtPb4ICvJwZbiMHduK2Z/dS7XVeWytSRMWorCeREREZlozuGWcy6W+BC4O9hm+1szMzONExcRmSPnXBAwJUNmKsXAxqLQjNcMjsRo6R6ktWeIlu7B8a1nYOz4icsdNHcNMjgSn1N/khzkZaWRnz1xK8hOIy84rszLYntpWPXERNYwM6Otd2iswPtokfdj9V20ByOykpMcW4pD3LmtmD1BofedZREy05KXuPciIiKyEszn28bkBEu/thcRWabSU5KpzMuiMm/mEWJmRs/gCP3DMQaH4wyOxBkciTE0MnrsR4uNHnf1D9PeO0Rb7xAdvUO09w5xuqmbjr5hOvqGJkyfdA42FGazsyzCzvIIu8r9l9aicPoiv3sRuZZGR2KdaeqeEGSdbemZMIo0Ky2ZTUUhnrtjPMjaURYhI1VBloiIiMzPfMKth4G5FeoSEZFlzTlHOCOVcEbqVb9WLG5E+4dp7x3kfEsvxxu6xlaR/NbhhrF2ReF0dpVH2FkWYXdFDs9Yl0dJJOOq7y8ii28kFudcSy/H6qMcr+/iWH0Xxxu6iPaPh1h5WalsLg7xwt2lbC4Oj61WWBbJ0JRmERERWVBzLigvM1NBeRGR6UX7hn3YFQRex+qjnG3uYSTu/y2qyM3kGevyxrbtpWFSVPxeZEkNDMc43hAEWPVRjtV3cbKxm6FgKnN6ShLbS8PsLM9hZ1mYrSU+yCoIaXSmiIiILJwFXS1xATpTamaN1/Sm15DCLRGRuRkciXGioZsnLnXw5KUOHr/UTlPXIACZqcnsq8oZC7v2VebqC7PIIorHjXMtPRys6eRgTSeHajs52dA9FkDnZKayqzwSbDnsLI+wsTBbIbSIiIgsuiUPt5xzKcCvAW8Bnm9mq/abicItEZGrM1q3ZzTseuJSB8cbuogFX67zslLZWBRiU1E2G4tCbCzMZlNxiOr8LFL1BVtk1syMpq7BsRDr4OVOjtRF6RkcASCcnsLeqhz2V+WytzKX3RU5lOdkzHYBIREREZEFtWThlnNuHz7Qeh1QgC8+b2a2aiuGKtwSEVl4fUMjHK6NcqQ2yvnWHs619HK+pYfWnqGxNilJjur8LDYWhdhQmEVJJIPiSAYl4XSKIxkUh9O1aqOsWQPDMc409XCisYuTDd2cbPRTC0dXK0xNduwoi7CvMpf9Vbnsq8plY2G2amOJiIjIsjFTuLXgP+U75/KB1+NDrX2jpxOadC30PUVEZHXLSkvh5o0F3LyxYML5aP8w51vGw67zLb2ca+nh4dMtDMXiT3ud7LTksaBrdJ+fnUZeVhp5WankjR5np5KbmUZaikaCycoSjxt1nf2cbOzmZIMPsE40dnGxtZdg8CMZqUlsKwnzvB0lbC8Ls68ql51arVBERERWsAUJt5wfn/5CfKD1EiCViYHWCPAD4HPANxbiniIiIjmZqVxXncd11XkTzpsZnX3DNHcP0tw9QHPX4KTjAQ7XdtLcNUj/cGza1w+lp5CXnUpeVhpZaclkpCaTkZJMRmqSP05NJj01KTjnz0cyUn1glp1GflYa+aE0stOSNZVLFlxrzyCnGrvHt6ZuzjR10zs0/v/0uoIstpeGefHecnaUhtleFqE6P4tkjcgSERGRVeSqwi3n3DbgzcAbgbLR08HegDPAPwP/aWYtV3MvERGR2XLO+VFY2WlsKw3P2HZgOEZH3xAdvcN+3zdER+8QHX3DtPcO0dk3RHvfMANDMdp7hxgYjjEwHA/2MQZG4mOrxk0nLTmJvOxU8rPTyQ/Csqr8LPZU5LCnIofKvEyFXzKtgeEYp5u6ORGMxBoNs9p6x6fl5mensa0kzCsPVLGtNOy3krCm4oqIiMiaMOefeJxzYeDV+FFaNyc+FezrgIrg+D/M7P9eVQ9FREQWUUZqMmU5mZTlZM77NeJxYygWp38oRteAD8U6+oZo6/H79t5h2nsHaQ8CtGP1XXz/WCPDsfEV6HZXRNgdhF17KnKozs9S4LXGmBmNXQOcaOjiREN3sO/iQsKUwszUZLaWhrlrRzHbSiNsLw2ztSRMUXjVrtUjIiIickWzDrecc8/GB1q/Dox+Axj9qbsX+BrwWeDH+GmIIiIia0JSkiMjyU9NzMtOY11B9hWvGRyJcaqxmyN1UY7WdXG0Lsonf3ZhLPAKZ6SwuzyHPZU57CqPsKs8hw2F2ZpOtkoMDMc429zD8SDAGh2V1dk3PNamMi+THWURXhRMKdwRTClUkXcRERGRiWYMt5xz6/DTDt8ErBs9Hezj+CDrs8BXzawv4boF76iIiMhqkp6SzN7KXPZW5o6dGxqJc7rJB14+9Iry6Z9fHCuOn5mazM7yCLvKI+wuz2FneYStJWEVvl/GzIyW7sEgxBofjXW+tZdYMBwrIzWJbaURXri7lB1lEXaURdhWGiaSkbrEvRcRERFZGa40cut8sE9Mq47hC8N/3szqF6VXIiIia1BaShK7K3LYXZHDa4Nzw7E4Z5t7OFoX5Vh9F8fqo3z1iVo++4tLAKQmO7aWhNlZFmFDUTYbC7PZUBhiXUGWVr+7xkanFR6qiXK4tpMjdVGO13dNqI1VnpPBjrIIL9g1GmSFWVegEXkiIiIiV+NK4ZbDF4Y34AvAP5jZwcXu1Fw45/KBLwHrgYvAq8ysY4p2dwMfAZKBj5vZB2e63jmXCnwcuB7/5/RZM/vrRX47IiIiE6QmJ42N5nllcC4eNy61900IvH5yqoX/eqJ27DrnoDwnkw2F2eNbEH5V5Wlq20Jo6xnkcF2Uw0GYdag2SmvPIAApST50vGtH8dh/v+2lYXKz0pa41yIiIiKrjzOz6Z90Lo4PtgCGge8BnwG+ZWbDs7ju/Wb2VwvX3Snv9SGg3cw+6Jx7N5BnZu+a1CYZOA08D6gFHgNea2bHp7veOfc64NfM7DXOuSzgOHCnmV2cqT8HDhywxx9/fMHfp4iIyJV0DwxzsbWP8609XGzt40JrDxdaeznf2kv3wHg5zOy0ZLYFNZxGRw9tL41oZb1JzIxo/zC1Hf3Ud/ZT1+n3Ne39HKmLUtfZD/ggcVNRiL2VOeyrzGVPZQ47yyIaOSciIiKygJxzT5jZgameu9JPsZ8CXgGEgTTgJcHW6Zz7MvA5M3tkITs7D/cCdwbHnwEeBN41qc2NwFkzOw/gnPticN3xGa43INs5l4IvoD8EdC3OWxAREbl64YxU9lT6IvSJzIz23iEutPZytrmHk43dHG/o4puH6vnCo5fH2q0ryGJH6Xjgta00vKpHefUOjtDYNUBjNNi6Bqjr7KcuIczqG4pNuCY9JYmKvEz2V+fypmeuY29lLrsrcggpGBQRERFZMjP+JGZmb3POvRN4Jb6w/O34qYp5wH3Afc65C4zX4Dq3uN2dUomZNQT9bXDOFU/RpgKoSXhcC9x0heu/gg++GoAs4A/NrH2qDjjn7sP/eVBdXX2Vb0dERGRhOecoCKVTEErnwPr8sfNmRl1n/4RC5ycauvjescaxNhmpSWwqCrG1JMzmYr/fUhyiKj9rWdeJiseN+mg/F1v7uNzeR2O0n4YgwBoNshJHs43Kz06jPDeDjUXZ3LalkIrcTL/l+X1+dpoWzhERERFZZq74a8ZgFcTPAJ9xzm0A3gL8BjCa4mwA3gu81zn3C+DzC91J59yPgNIpnvrz2b7EFOemn4/p3QjEgHJ8mPdT59yPRkd/TXghs/uB+8FPS5xln0RERJaUc47KvCwq87J43s6SsfO9gyOcbOzmbHM3p5t6ONPcwy/Pt/G1p+rG2qSnjIZeIcpyMynITqMonE5BdjqF4TQKstPJz05b1ADMzGjqGuRCay8X23q5GEzBvNjay6X2PoZG4mNtkxwUhdMpjfjg6pmbCijNyaQ0J52SSAZlOZmURjLITNNUQhEREZGVZk5j6M3sAuNB1l34oOtl+Gl7ALcE26j1zrl0Mxu8mk6a2XOne8451+ScKwtGXZUBzVM0qwWqEh5XAqMrPU53/euA7wW1xZqdcz8HDjC+gqSIiMiqlJ2ewjPW5fGMdXkTzncPDHOmuYezTT2cCYKvxy520NzdwHDs6b/bcQ7ys9IoCKVRGEonlJ5CRmoy6SlJU+9Tk8hISSZuRs/giN8GRugdGqF7YITe0XODMXoGh2nrGZowbTAtOYnqgiw2FGbz7O3FrC/IZn1hFusLsikOp5OSnLTof3YiIiIicu3Nu0CEmT0APOCci+CDoDfjRzvB+KioNwG/7pz7Cn7a4oPz7+q0vhnc54PB/htTtHkM2BKMPKsDXhP0eabrLwPPcc59Hj8t8Wbgw4vQfxERkRUhnJHK9dV5XF89MfQyM7r6R2jtHaS1e5C23iFaewZp7fH7tuC4vbePwZE4A8Oxsf3AcIz4DGOes9KSCaWnEEpPITvYV+SmEc4Ik5uVyobCbNYX+NUgy3Mzl/VUSRERERFZHDOuljjnF3NuB/BW4PVMnEY4epM64Atm9p4FvGcB8GX8NMnLwCvNrN05Vw583MzuCdrdgw+nkoFPmtkHrnB9CF9Qfyd+WuOnzOxvr9QfrZYoIiIyN8Ox+ITQywGhjBSy01IUVomIiIgIMPNqiQsabiXcMBl4IX7a4ouB1ISnzcxWbUELhVsiIiIiIiIiIgtrpnBrUYpPmFnMzL5lZi/HF2T/f4DDi3EvERERERERERFZuxa9sqqZtZnZh81sP74g+8cW+54iIiIiIiIiIrI2zLug/HyY2ZPAk9fyniIiIiIiIiIisnppTWwREREREREREVmxFG6JiIiIiIiIiMiKtSirJa5lzrkW4NJS92OBFAKtS90JkRVAnxWR2dFnRWR29FkRmR19VkRmZ7V8VtaZWdFUTyjckmk55x6fbplNERmnz4rI7OizIjI7+qyIzI4+KyKzsxY+K5qWKCIiIiIiIiIiK5bCLRERERERERERWbEUbslM7l/qDoisEPqsiMyOPisis6PPisjs6LMiMjur/rOimlsiIiIiIiIiIrJiaeSWiIiIiIiIiIisWAq35Gmcc3c75045584659691P0RWS6cc1XOuZ8450445445534/OJ/vnPuhc+5MsM9b6r6KLAfOuWTn3FPOuW8Fj/VZEZnEOZfrnPuKc+5k8O/LLfqsiDydc+4Pg5+/jjrn/tM5l6HPigg45z7pnGt2zh1NODftZ8M5957gu/4p59wLlqbXC0/hlkzgnEsGPga8ENgJvNY5t3NpeyWybIwAf2RmO4Cbgd8JPh/vBh4wsy3AA8FjEYHfB04kPNZnReTpPgJ8z8y2A/vwnxl9VkQSOOcqgN8DDpjZbiAZeA36rIgAfBq4e9K5KT8bwXeX1wC7gmv+OcgAVjyFWzLZjcBZMztvZkPAF4F7l7hPIsuCmTWY2ZPBcTf+C0gF/jPymaDZZ4CXLkkHRZYR51wl8CLg4wmn9VkRSeCciwC3A58AMLMhM+tEnxWRqaQAmc65FCALqEefFRHM7GGgfdLp6T4b9wJfNLNBM7sAnMVnACuewi2ZrAKoSXhcG5wTkQTOufXAdcCjQImZNYAPwIDiJeyayHLxYeBPgXjCOX1WRCbaCLQAnwqm8H7cOZeNPisiE5hZHfB3wGWgAYia2Q/QZ0VkOtN9Nlbt932FWzKZm+KcltQUSeCcCwFfBf7AzLqWuj8iy41z7sVAs5k9sdR9EVnmUoDrgX8xs+uAXjStSuRpgnpB9wIbgHIg2zn3hqXtlciKtGq/7yvckslqgaqEx5X4Ib8iAjjnUvHB1hfM7L+D003OubLg+TKgean6J7JM3Ar8mnPuIn56+3Occ59HnxWRyWqBWjN7NHj8FXzYpc+KyETPBS6YWYuZDQP/DTwTfVZEpjPdZ2PVft9XuCWTPQZscc5tcM6l4YvNfXOJ+ySyLDjnHL4uygkz+4eEp74JvCk4fhPwjWvdN5HlxMzeY2aVZrYe/+/Ij83sDeizIjKBmTUCNc65bcGpu4Dj6LMiMtll4GbnXFbw89hd+Nqn+qyITG26z8Y3gdc459KdcxuALcCvlqB/C86ZrYoRaLKAnHP34GulJAOfNLMPLG2PRJYH59xtwE+BI4zXEfozfN2tLwPV+B++Xmlmk4s6iqxJzrk7gT82sxc75wrQZ0VkAufcfvzCC2nAeeAt+F9A67MiksA5937g1fjVq58C3g6E0GdF1jjn3H8CdwKFQBPwl8DXmeaz4Zz7c+Ct+M/SH5jZd699rxeewi0REREREREREVmxNC1RRERERERERERWLIVbIiIiIiIiIiKyYincEhERERERERGRFUvhloiIiIiIiIiIrFgKt0REREREREREZMVSuCUiIiKyBjjn1jvnLNg+vUCveTF4vYsL8XoiIiIi85Gy1B0QERERWc6cczbD0z1AE/Ak8N/AV81s+Jp0bBE4594MrAcws/ctZV9EREREZkvhloiIiMj8hYJtE/BK4Ihz7hVmdnppuzVvbwbuCI7ft3TdEBEREZk9ZzbTLyNFRERE1rZJI7deNunpPOCZwOuBzODcZeA6M2u/Bt1bUM65BwnCLTNzS9sbERERkdlRuCUiIiIyg8Rwa7rAxzm3C3gQKAxOfcjM3rX4vVtYCrdERERkJVJBeREREZGrZGbHgD9LOPWKpeqLiIiIyFqjcEtERERkYXw74Xijcy5rcgPn3A7n3Eecc0edc1HnXL9z7pJz7svOuclTHqfknLveOfevzrkjzrku59ywc67ZOXfcOfc/zrl3Ouc2THHdtKslOuceDEao3ZFwzqbY3jfpulmvluice75z7nPOufPOuT7nXLdz7mTwXp5xhWuf1nfnXKFz7n3Bn0N3sD3pnHvPVH/2IiIisnqpoLyIiIjIwmiZ9DgX6Bt94Jx7P/DnQPKkdtXB9krn3EPAy82sbaobBOHSe4HJUwaLgm0H8GLgLuCl83gPC845FwL+A3jJFE9vC7b7nHP/BPyhmcVn8ZoHgK8DFZOeui7YXuWcu2sl1j0TERGRuVO4JSIiIrIwiiY97ho9cM79NfDu4GEM+CLwY6Af2AO8FSjBj5z6sXPuZjPrT3wx59y9wF8GD/uB/wR+CbQDGUAlcAB43jz6/r/x9cL+X2BXcG6qkWQn5/Kizrlk4LvAbcGpTuCTwJP4n0NvA34DSAN+D1+U/74rvGwVfpRcPvAF4CdAD7AT+B2gANgPfDh4bREREVnlVFBeREREZAazKSgftPtN4P7g4UUz2xCcvwX4OX60VS9wj5k9POnafOD7+HAK4O/M7E8mtfkW8CJ8OHa7mT0yTT8ygL1m9qtJ59cDF4KHnzGzN09x7YPMoaB8MB1xHXDJzNZP8fy7gA8GD08BzzGz+kltrgN+hA+rAF5iZt+aoe/gQ7K7zezRSe024IOzXPyfU/Xk+4mIiMjqo5pbIiIiIlfJObcd+EDCqa8kHP8J49MI/2RysAUQTJ97BePTGH/bOZc7qdnmYH9sumAreK2BycHWUnDOpQF/GDwcAV45VdBkZk8Bv5Vw6t2T20zh9yYHW8FrXQA+FjxMxk/PFBERkVVO4ZaIiIjILDnnXjppe7Nz7n7gCcanJdYDHwrapwP3BOfbgE9M99pmdgk/1RAgG3j+pCajwVelcy7n6t/NonsmfqolwHfN7Mh0Dc3sK8DZ4OGtzrniGV63BV/Dazo/TjjeOZuOioiIyMqmcEtERERk9r42afsU8JvA6Op8J4Dnm9locfl9QHpw/KCZDV3h9X+QcHzTpOd+GOzzgYecc691zkXm/haumRsTjn8wbatxP0w4nvzeEz1uZrEZnq9LOM6bxX1FRERkhVNBeREREZH56wWagafwYdd/mdlgwvNlCcenZ/F6iW3KJj33QfxKiDvxodl/ADHn3EF8Ta+fAN+fXIh+CS3ke0/UeoXXSfzzz5jFfUVERGSFU7glIiIiMkuzKbI+STjhuHcW7XumuRYz63DO3Qy8C3g7fspfMvCMYPs9oNs592Hg/53FKLHFtmDvfZL4/LojIiIiq5WmJYqIiIgsnu6E4+xZtA9Ncy0AZtZtZv8bKAeuB94JfInx0Uxh4C+Abzrn5hrELbQFfe8iIiIi01G4JSIiIrJ4GhKOt8yifWKbp60sOMrM4mb2lJl91Mxegx/F9TKgPWjyAuBFc+3sAluU9y4iIiIymcItERERkcVziPEaUHc651Kv0D5xhcRfzfYmQdj1deC9Cadvm+31Ccam/C3AyK/E/j9vFu0T28z6vYuIiIgo3BIRERFZJEFx+W8HDwuBN0/X1jlXBbw2eNjL7FYYnOxiwvF8aqsm1r2azVTCmTwCNAbHL3LO7ZyuoXPu1xkfufUzM2u+ynuLiIjIGqJwS0RERGRx/S3jI6L+3jl36+QGzrk84CuMB0r/Ymadk9rc75zbPd1NnHMpwG8mnDo0j75eSDi+fh7XjwkK2v9j8DAF+C/n3NNWQXTO7QX+LeHUB6/mviIiIrL2aLVEERERkUVkZr90zv0N8B58wfeHnHP/CfwY6Ad2M776IcBhJk4vHPWbwG86544BPwGO4mtsZQMbgdcwPvrpND4sm6sH8KsuAnzCOfePwCUgFpw7a2Zn5/B6fw+8BD9FcidwzDn3SeBJ/M+htwJvAtKD9v9uZt+e6oVEREREpqNwS0RERGSRmdmfOedGgD8DkoE3BNtkDwEvN7P+qV4GcMCuYJvOYeDeaV7jSr4N/AwfRm0GPjbp+fcD75vti5lZzDn3QuA/gRcDecAfTdU0uNfvz73LIiIistZpWqKIiIjINWBm7wX2Av8EHAe68cXma4Gv4kOtO82sbZqXKMXX5Pp3/MinDvyIqn58ra2vAa8Hrjezi/PsYwxf2P3dwC8S7jFvZtZjZi8B7gb+Az8SbABfV+w0cD9wg5m908zi07+SiIiIyNScmS11H0REREREREREROZFI7dERERERERERGTFUrglIiIiIiIiIiIrlsItERERERERERFZsRRuiYiIiIiIiIjIiqVwS0REREREREREViyFWyIiIiIiIiIismIp3BIRERERERERkRVL4ZaIiIiIiIiIiKxYCrdERERERERERGTFUrglIiIiIiIiIiIrlsItERERERERERFZsRRuiYiIiIiIiIjIivX/A8BCQwCo4mIVAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.figure(figsize=(20,3))\n", "plt.plot(sol[500])\n", "plt.plot(sol[10000])\n", "plt.xlabel('Position', fontsize=30)\n", "plt.ylabel('Amplitude', fontsize=30)\n", "plt.savefig('../output/sample.png')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Make an animation of the string: (1/dt is fps)" ] }, { "cell_type": "code", "execution_count": 117, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "50000" ] }, "execution_count": 117, "metadata": {}, "output_type": "execute_result" } ], "source": [ "len(sol[::10,:])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So 200000 fps. Lets only index every 10 frames which gives 20000 fps. If we choose a gif fps of 20, this means our animation is moving 1000x slower than real life" ] }, { "cell_type": "code", "execution_count": 118, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAD8CAYAAACPWyg8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAoCUlEQVR4nO3deXxV1b338c8vCWEeJUBIwiRhlEkiaFUcUUAtzoLjpVpqq61tn94Wb+vT3tvh+vS2vdXWqjhUaFUcagsiFoU6V4QgyIyEOSEkQSAMgYQkv+ePs2nTmOFsciDT9/16ndc5e++19l6LId/svfZex9wdERGRMOLquwEiItL4KDxERCQ0hYeIiISm8BARkdAUHiIiEprCQ0REQotJeJjZBDPbaGZZZjajiu1mZg8H21eZ2ZkVtj1tZvlmtqZSnS5m9qaZbQreO1fYdn+wr41mdnks+iAiItGrc3iYWTzwCDARGAJMNbMhlYpNBNKD13Tg0QrbngEmVLHrGcBid08HFgfLBPueAgwN6v0uaIOIiJwisTjzGANkufsWdy8B5gCTK5WZDMz2iCVAJzNLBnD3d4G9Vex3MjAr+DwLuLrC+jnuXuzuW4GsoA0iInKKJMRgHynAzgrL2cDYKMqkALk17Le7u+cCuHuumXWrsK8lVezrc8xsOpEzHdq2bTt60KBBNfdERET+xfLly/e4e1Ll9bEID6tiXeU5T6IpE8vjRVa6zwRmAmRkZHhmZuYJHlJEpHkys+1VrY/FZatsIK3Cciqw6wTKVJZ3/NJW8J5fh32JiEgMxSI8lgHpZtbXzBKJDGbPq1RmHnB7cNfV2UDh8UtSNZgH3BF8vgOYW2H9FDNraWZ9iQzCL41BP0REJEp1vmzl7qVmdi+wEIgHnnb3tWZ2d7D9MWABMInI4HYRMO14fTN7HrgQ6Gpm2cAP3f0p4EHgRTO7E9gB3BDsb62ZvQisA0qBe9y9rK79EBGR6FlzmZJdYx4iIuGZ2XJ3z6i8Xk+Yi4hIaAoPEREJTeEhIiKhKTxERCQ0hYeIiISm8BARkdAUHiIiEprCQ0REQlN4iIhIaAoPEREJTeEhIiKhKTxERCQ0hYeIiISm8BARkdAUHiIiEprCQ0REQlN4iIhIaAoPEREJLSbhYWYTzGyjmWWZ2YwqtpuZPRxsX2VmZ9ZW18xeMLOVwWubma0M1vcxsyMVtj0Wiz6IiEj0Euq6AzOLBx4BxgPZwDIzm+fu6yoUmwikB6+xwKPA2JrquvtNFY7xS6Cwwv42u/vIurZdREROTCzOPMYAWe6+xd1LgDnA5EplJgOzPWIJ0MnMkqOpa2YG3Ag8H4O2iohIDMQiPFKAnRWWs4N10ZSJpu75QJ67b6qwrq+ZrTCzd8zs/Lo0XkREwqvzZSvAqljnUZaJpu5U/vWsIxfo5e6fmdlo4C9mNtTdD3yuYWbTgekAvXr1qqb5IiISVizOPLKBtArLqcCuKMvUWNfMEoBrgReOr3P3Ynf/LPi8HNgMDKiqYe4+090z3D0jKSkpZLdERKQ6sQiPZUC6mfU1s0RgCjCvUpl5wO3BXVdnA4XunhtF3UuBDe6efXyFmSUFA+2YWT8ig/BbYtAPERGJUp0vW7l7qZndCywE4oGn3X2tmd0dbH8MWABMArKAImBaTXUr7H4Knx8oHwf8l5mVAmXA3e6+t679EBGR6Jl75SGGpikjI8MzMzPruxkiIo2KmS1394zK6/WEuYiIhKbwEBGR0BQeIiISmsJDRERCU3iIiEhoCg8REQlN4SEiIqEpPEREJDSFh4iIhKbwEBGR0BQeIiISmsJDRERCU3iIiEhoCg8REQlN4SEiIqEpPEREJDSFh4iIhKbwEBGR0BQeIiISWkzCw8wmmNlGM8sysxlVbDczezjYvsrMzqytrpn9yMxyzGxl8JpUYdv9QfmNZnZ5LPogIiLRS6jrDswsHngEGA9kA8vMbJ67r6tQbCKQHrzGAo8CY6Oo+7/u/otKxxsCTAGGAj2BRWY2wN3L6toXERGJTizOPMYAWe6+xd1LgDnA5EplJgOzPWIJ0MnMkqOsW9lkYI67F7v7ViAr2I+IiJwisQiPFGBnheXsYF00ZWqre29wmetpM+sc4ngAmNl0M8s0s8yCgoJo+yMiIrWIRXhYFes8yjI11X0UOB0YCeQCvwxxvMhK95nunuHuGUlJSVUVERGRE1DnMQ8iv/mnVVhOBXZFWSaxurrunnd8pZk9AcwPcTwRETmJYnHmsQxIN7O+ZpZIZDB7XqUy84Dbg7uuzgYK3T23prrBmMhx1wBrKuxripm1NLO+RAbhl8agHyIiEqU6n3m4e6mZ3QssBOKBp919rZndHWx/DFgATCIyuF0ETKupbrDrn5vZSCKXpLYBXwnqrDWzF4F1QClwj+60EhE5tcy9yuGCJicjI8MzMzPruxkiIo2KmS1394zK6/WEuYiIhKbwEBGR0BQeIiISmsJDRERCU3iIiEhoCg8REQlN4SEiIqEpPEREJDSFh4iIhKbwEBGR0BQeIiISmsJDRERCU3iIiEhoCg8REQlN4SEiIqEpPEREJDSFh4iIhKbwEBGR0GISHmY2wcw2mlmWmc2oYruZ2cPB9lVmdmZtdc3sf8xsQ1D+z2bWKVjfx8yOmNnK4PVYLPogIiLRq3N4mFk88AgwERgCTDWzIZWKTQTSg9d04NEo6r4JnOHuw4FPgfsr7G+zu48MXnfXtQ8iIhJOLM48xgBZ7r7F3UuAOcDkSmUmA7M9YgnQycySa6rr7m+4e2lQfwmQGoO2iohIDMQiPFKAnRWWs4N10ZSJpi7Al4DXKyz3NbMVZvaOmZ1fXcPMbLqZZZpZZkFBQe09ERGRqMQiPKyKdR5lmVrrmtn3gVLg2WBVLtDL3UcB3waeM7MOVTXM3We6e4a7ZyQlJdXQBRERCSMhBvvIBtIqLKcCu6Isk1hTXTO7A7gSuMTdHcDdi4Hi4PNyM9sMDAAyY9AXERGJQizOPJYB6WbW18wSgSnAvEpl5gG3B3ddnQ0UuntuTXXNbALwPeCL7l50fEdmlhQMtGNm/YgMwm+JQT9ERCRKdT7zcPdSM7sXWAjEA0+7+1ozuzvY/hiwAJgEZAFFwLSa6ga7/i3QEnjTzACWBHdWjQP+y8xKgTLgbnffW9d+iIhI9Cy4GtTkZWRkeGamrmyJiIRhZsvdPaPyej1hLiIioSk8REQkNIWHiIiEpvAQEZHQFB4iIhKawkNEREJTeIiISGgKDxERCU3hISIioSk8REQkNIWHiIiEpvAQEZHQFB4iIhKawkNEREJTeIiISGgKDxERCU3hISIioSk8REQktJiEh5lNMLONZpZlZjOq2G5m9nCwfZWZnVlbXTPrYmZvmtmm4L1zhW33B+U3mtnlseiDiIhEr87hYWbxwCPARGAIMNXMhlQqNhFID17TgUejqDsDWOzu6cDiYJlg+xRgKDAB+F2wHxEROUUSYrCPMUCWu28BMLM5wGRgXYUyk4HZ7u7AEjPrZGbJQJ8a6k4GLgzqzwLeBr4XrJ/j7sXAVjPLCtrwYQz68jn/+epadhce/cfyRQO7ceNZaSfjUCIiJ8zd2ZR/iDfX5bF2VyHu/9z20JRRJCbEdpQiFuGRAuyssJwNjI2iTEotdbu7ey6Au+eaWbcK+1pSxb4+x8ymEznToVevXlF251/t3HuEHXsPA3C4uIzX1+wm78BR7r24P2Z2QvsUEYmVQ8Wl/OHD7bywbAfbPisCoF/XtiTE//Pnk+PVVT9hsQiPqn6CVm5pdWWiqXsix4usdJ8JzATIyMg4oT+9J+/I+Mfn0rJyvvunVfzyzU8pOlbGdy8fqAARkXpxqLiUZz7YypPvb2V/0TG+cPppfHlcPy4d3J3uHVqd9OPHIjyygYrXcVKBXVGWSayhbp6ZJQdnHclAfojjnRQJ8XH84voRtG4Rz6Nvb6aktJwfXDFYASIip4y78+qqXH4yfx35B4u5eFA3vnFJOiPTOp3SdsQiPJYB6WbWF8ghMph9c6Uy84B7gzGNsUBhEAoFNdSdB9wBPBi8z62w/jkz+xXQk8gg/NIY9CMqcXHGT64+gxbxcTz1/lZ6dGjFl8f1O1WHF5FmLCv/ID+at473s/YwLKUjj982mlG9Otde8SSoc3i4e6mZ3QssBOKBp919rZndHWx/DFgATAKygCJgWk11g10/CLxoZncCO4AbgjprzexFIoPqpcA97l5W136EYWb83yuHUHCwmJ8uWE9yp1ZcObznqWyCiDQj2/Yc5qHFm5i7Moe2LRP48eSh3Dy2N/Fx9XfVw9xjP5DSEGVkZHhmZmZM93n0WBm3PfURn+ws5I93jWVM3y4x3b+INF/uzvLt+/jjku28uiqXFvHGHef0Yfq4fpzWruUpa4eZLXf3jM+tV3jUzf6iEq599O98dqiEl+4+hwHd28f8GCLSfJSWlfOnj7N5+v1tbMw7SPuWCdyQkcbdF/ajW/uTPxBemcLjJIUHwM69RVz76N9JiDNe+doXSO7Y+qQcR0SaLndn8fp8HvzrBrLyDzG0ZwduO7s3V43oSduWsRiePjHVhYfmtoqBtC5teGbaWRw8WsodTy+lsOhYfTdJRBqRjbsPMvWJJdw1O5PycuexW0cz/+vnMWVMr3oNjpooPGJkaM+OzLxtNFv3HGbaM0s5cFQBIiI1KzxyjP98dS2THn6PDbsP8uPJQ1n4rXFMOKNHg38EoGFGWiP1hf5deXjKKL4xZwU3P7GEWdPGnNKBLRFp+I6UlPHWxnxeW5XL3zbkc7S0jJvH9OI7lw2kc9vE+m5e1BQeMTZxWDIzE+O5+w/LufHxD3n2rrPp0fHUD3KJSMNx/M6plzKzeW11LoeKS+naLpHrRqcw5axenJHSsb6bGJoGzE+Sj7Z8xp2zMunYugVP3pHB4OQOp+zYItJwvL9pDz95bR0bdh+kTWI8k4Ylc82oFMb27UJCfMMfOdDdVqc4PADW5BRy56xlHDpaykNTRnHpkO6n9PgiUn+27TnMTxes5811efTq0oZ7LjqdK4b3pF0DHQCvjsKjHsIDIO/AUb48O5PVOYV8f9Jg7jpfU5mINGWHi0v57VtZPPXeVlrEG/dc3J8vnduXVi0a59cOVRcejSsCG6HuHVrxwvRz+M5Ln/CT19ZT7s70cafXd7NEJMbcnXmf7OJnC9aTd6CYa89MYcaEQXQ7BTPc1geFxynQOjGeh6eOwgx+tmADrRMTuO3s3vXdLBGJkY937OPH89exYsd+hqV05He3jGZ07/qZsPBUUXicIvFxxv/eNJKjx8p44C9raNMinutGp9Z3s0SkDnbuLeIXb2xk7spdJLVvyc+vH851Z6bW64SFp4rC4xRqER/Hb28+k7tmZfLvL39C25YJTDijR303S0RC2l9Uwm//lsXsD7djBvde1J+vXnh6g30a/GRoPj1tIFq1iOfx20Zz61Mf8Y3nV/D0v53Feeld67tZIhKFvYdL+P0HW3nm79s4VFzK9Wem8u3LBjTL+ex0t1U9KSw6xk0zP2TH3iL+eNdYzqynL3QRkdptLjjE8x/t4LmlOygqKWPC0B58c3w6g3o0/ee3dKtuAwsPgPyDR7nxsQ8jv81MG9PkB9hEGpPi0jJe+TiHlzJ38vGO/cTHGVcNT+ZrF/VvVl+9oPBogOEBkL2viFuf/Ii8A8XMvH0056cn1XeTRJq10rJyXlmRw0OLNpGz/wgDurfjhtFpTB7Vs16+T6O+KTwaaHhA5Azk9qeWsqXgMA9PHcmEM5Lru0kizU55ufPa6lx+vehTNhccZkRqR/798kGc2/+0Bj/D7cl0Ur7Pw8y6mNmbZrYpeK/yuouZTTCzjWaWZWYzaqtvZuPNbLmZrQ7eL65Q5+1gXyuDV7e69KEh6NY+8iDhGSkd+NqzHzN/1a76bpJIs+HuLFidy4SH3uXrz68gzozHbh3NX+45l/PSuzbr4KhJXWflmgEsdvd0YHGw/C/MLB54BJgIDAGmmtmQWurvAa5y92HAHcAfKu32FncfGbzy69iHBqFjmxb84c6xjO7dmfvmrOSva3Lru0kiTV7O/iPc/vRSvvbsx5SVO7+ZOoqF32wc36dR3+oaHpOBWcHnWcDVVZQZA2S5+xZ3LwHmBPWqre/uK9z9+K/fa4FWZtbkvxijbcsEfj9tDCNSO3Lvcyt4c11efTdJpElyd+Ys3cHl//suy7fv48eTh/LGty7gqhE9iWsGD/jFQl3Do7u75wIE71VdQkoBdlZYzg7WRVv/OmCFuxdXWPf74JLVA1bDrwdmNt3MMs0ss6CgIPpe1aN2LRN45ktjGNqzA197drkuYYnE2L7DJUz/w3JmvLKaYSkdWfjNcdx2Tp9m8VR4LNX6kKCZLQKqegz6+1Eeo6q/kahG6c1sKPD/gMsqrL7F3XPMrD3wJ+A2YHZV9d19JjATIgPmUba33nVo1YLZd47lrlnL+PrzK9h3uITbzulT380SafSWbt3LfXNWsOdQMT+4YjBfOrevzjROUK3h4e6XVrfNzPLMLNndc80sGahq/CEbSKuwnAoc/3W62vpmlgr8Gbjd3TdXaE9O8H7QzJ4jclmsyvBozDq2joyB3Pvcxzwwdy0Fh0r41qXpug4rcgLyDx7l4cWbeO6jHaR1acMrXz2XYamN79v7GpK6XraaR2RAm+B9bhVllgHpZtbXzBKBKUG9auubWSfgNeB+d//g+I7MLMHMugafWwBXAmvq2IcGq1WLeB67dTQ3jE7l4cWbeGDuGsrKG80JlEi92114lF+9sZELfv42c5bu5JaxvZn/9fMUHDFQ17mtHgReNLM7gR3ADQBm1hN40t0nuXupmd0LLATigafdfW1N9YF7gf7AA2b2QLDuMuAwsDAIjnhgEfBEHfvQoCXEx/Hz64fTpV0ij7+zhX1Fx/jVjSNomdA4v1hG5GTLO3CUlzJ38sa6PFZlFwJwxfBkvnPZQPp2bVvPrWs69JBgIzLz3c38bMEGzuvflcduG93ovs5S5GT67FAxj72zmdkfbqe4tJyRaZ0YP6Q7lw/tQf9u7eq7eY2WvkmwCZg+7nS6tG3J9/60ilueWMIz08bQuW1ifTdLpF6VlpXz1PtbeXjxJo4cK+PqUSncd0k6vU/TWcbJpPBoZK4fnUrH1i2457mPueHxD/nDnWOa5XTQIgAbdh/guy+vYlV2IZcO7saMiYPo3635TFpYn+o6YC71YPyQ7sz+0hjyCo9y/aMfsiansL6bJHJKrc4u5D/+vJqrfvM+OfuO8NubR/HE7RkKjlNIYx6N2JqcQu6ctYx9Rcd44Moh3Dq2l27llSbrUHEpc1fm8PzSHazJOUDLhDiuGZXCv18+kNPaNfkJKOqNZtVtguEBkUHCb7/4Ce98WsAVw5L5wZWDdRlLmpRP8w7y+w+2MW9lDodLyhic3IGpY9KYPDKFjq1b1HfzmjwNmDdRp7Vrye//7Swef3cLv3hjIwvX7uaLI3syfVy/ZvEtZ9J0bSk4xK8XbeLVVbtomRDHVcN7cvPYXoxM66Qz7AZAZx5NyM69RTz1/lZeWLaTI8fKGNqzAxOG9mDCGT1Ib0bffCaN2+rsQp56fwvzPtlFy4R47vhCH74yrp/uLKwnumzVDMLjuP1FJby8PJvX1+xm+fZ9AJyf3pXvXj5IT9ZKg1Re7ixan8eT721l6ba9tE2MZ+qYXnzlgtNJaq/xjPqk8GhG4VFR3oGj/HlFDo+/s5l9RceYeEYPvj1+gM5EpEEoKS1n7socHn93C1n5h0jp1Jpp5/bhxrPS6NBK4xkNgcKjmYbHcQePHuOp97fy5HtbOVxSyuQRPbnv0gGarkHqzdsb8/nRvLVs+6yIQT3a89ULT+eKYckkxOsJgoZE4dHMw+O4vYdLePzdzcz6+zaOlTkXDezG1aN6cung7rRqofmy5OTLLTzCj+evY8Hq3fTr2pYfXDmYiwZ20yB4A6XwUHj8i/yDR3nq/a38+eMc8g8W065lAhcMTOLCAUlcMDCJbu1b1XcTpYlxd+Ys28lPX1vPsbJyvnFJOned31eTfDZwCg+FR5XKyp0lWz5j3spd/G1jPgUHI1/YmNG7M9ecmcKVw3rSsY2uPUvd5Ow/wow/reK9TXs4u18Xfn7dCHqd1qa+myVRUHgoPGpVXu6syz3AWxvymfvJLrLyD5EYH8f4Id256aw0zuvfVd+6JqGUlzvPfrSdB1/fgAP3TxzELWN7699RI6LwUHiE4u6syTnAKyuy+cuKHPYVHSOlU2uuGZXCFcOTGdSjva5RS42y8g9y/yurWbZtH+f178p/XzuMtC4622hsFB4KjxNWXFrGm+vyeGHZTj7I2kO5Q7+ubblyeDJfHNlTk9HJP5SVO29tyOcPS7bzzqcFdGzdgh9cMZjrR6fql41GSuGh8IiJPYeKWbh2N6+tymXJls8odxic3IHhKR1xnHKPXKo4Vu6UlpXTqU0Lhqd2YnhqRwZ2b6/bMJuYsnJn6da9LNu2l1XZ+1mxYz+fHS6he4eWTB3Ti1vP7k1XTVrYqCk8FB4xl3/gKPNX5fLqql3s2n+EODMMiIszWsTHkRBn5B8spvDIMQA6tErg4kHduGxoDy4YkERbfRNio7W54BAvL49c0swtPIoZnJ7UjuGpHbl0cHfGD+lOC/2i0CSclPAwsy7AC0AfYBtwo7vvq6LcBOAhIt87/qS7P1hTfTPrA6wHNga7WOLudwd1RgPPAK2BBcB9HkUnFB71w93Z/lkRn2Tv5/1Ne1i0Po99RcdIjI9jTN8uXDSoGxcOTKJf17a6rNEIrNt1gN/8bROvr9lNfJxxwYAkrj0zhXEDkvREeBN1ssLj58Bed3/QzGYAnd39e5XKxAOfAuOBbGAZMNXd11VXPwiP+e5+RhXHXArcBywhEh4Pu/vrtbVV4dEwlJaVs2zbPt7amM9bG/LZlH8IgLQurblgQBLjh/RgXHpXBUkDs3z7Xh57ZwtvrsujfcsEpp3bh1vP6a3ngZqBkxUeG4EL3T3XzJKBt919YKUy5wA/cvfLg+X7Adz9v6urX114BGXecvdBwfLUoP5XamurwqNh2rm3iHc+LeDtjQX8ffMeikrKGJLcga9f3J/Lh/bQLZ31yN1ZvD6fx97ZTOb2fXRs3YJp5/Zh2rl99T0azcjJ+j6P7u6eCxAEQLcqyqQAOyssZwNjo6jf18xWAAeAH7j7e8G+sivtK6W6xpnZdGA6QK9evUJ1TE6NtC5tuPXs3tx6dm9KSst59ZNdPPJWFl999mMGdm/Pdy4fyKWDNXXFqZa5bS8/W7Cej3fsJ6VTa3541RBuzEjTOJX8Q63/EsxsEdCjik3fj/IYVf2vr+10Jxfo5e6fBWMcfzGzoWH35e4zgZkQOfOIsr1STxIT4rhudCpXj0ph/qpdPLRoE1+enUlG78585/KBjO3bRSFyku0uPMoP561h4do8urVvyYPXDuP60am6S04+p9bwcPdLq9tmZnlmllzhslN+FcWygbQKy6nAruBzlfXdvRgoDj4vN7PNwIBgX6nV7EuaiPg4Y/LIFCYNS+alzGx+vehTpsxcwoDu7bjprF5cMyqFLvpioJhyd+au3MX/nbuGkrJyvj1+AHed35c2iTrTkKrV9V/GPOAO4MHgfW4VZZYB6WbWF8gBpgA311TfzJKIDKSXmVk/IB3Y4u57zeygmZ0NfATcDvymjn2QBqpFfBw3j42ExZ9X5PBC5k5+PH8d/71gPRl9OnPxoG5cPKg7/bu1q++mNmp5B47yo3lreX3Nbs7s1Ylf3jhSU/VLreo6YH4a8CLQC9gB3BD8gO9J5JbcSUG5ScCvidyq+7S7/7SW+tcB/wWUAmXAD9391aBOBv+8Vfd14Ou6Vbf52LD7AHNX7uKtDfls2H0QiDykOHlkT64a0ZOUTq3ruYWNR0lpOc/8fSsPLdrEsTLnm+PT+cq404nXTQpSgR4SVHg0OTn7j/DG2t3M+2QXK3bsB2Bs3y5cPSqFSWckazbgapSXO39du5tfvrGRzQWHuWRQNx64cgh9dLYhVVB4KDyatO2fHWbeyl38eWUOWwoOkxgfx8WDunH1qBQuGpSk74wgMpXIm+vy+PWiT9mw+yD9ktry/UmDuWRw9/pumjRgCg+FR7NQcTbgVz/JZc+hYjq2bsE9F53Ov32hL4kJze+uoZ17i3hpeTYvZ+5kV+FR+nZtyzcu6c8XR6ToEpXUSuGh8Gh2SsvK+WDzZzzzwVbe2lhAv65teeCqIVw0sKrHkZqWvYdLWLA6l3mf7GLZtr0AnNe/KzedlcaEoT10661ETeGh8GjW3tqQz4/nr2PLnsOcn96V714+iGGpHeu7WTHj7mzMO8hbGwp4e2M+y7fvo7TcOT2pLVePTOHa0am6mUBOiMJD4dHslZSWM/vDbTzyVhb7io5xxfBkvnrB6ZyR0jhDpLi0jA+y9rB4fWSesF2FR4HI3WcXDkziyuHJDEnuoAcrpU4UHgoPCRw8eown3t3Ck+9vpaikjBGpHbnl7N58cURPWrVo2APr7s7aXQd4KXMncz/Zxf6iY7RNjOe89K5cPKgbFw7sRvcOmqxQYkfhofCQSgqPHOOVj7N59qMdZOUfomfHVnzjknSuG53a4L6LoqzceWPtbma+t4UVO/aTmBDHZUO6c93oVL5w+mm6m0xOGoWHwkOq4e68n7WHX7zxKZ/s3E+f09pw2zl9uGJYMj061u9v8UdKynh5+U6eeG8rO/YW0atLG750bh+uGZWq51jklFB4KDykFu7OovX5PLx4E6tzCjGDs3p34YKBSYzu3ZkRqZ1onXhqfsPfc6iYZ5fsYNaH29h7uIRRvTrxlXH9GD+kh26vlVPqZE3JLtJkmBnjh0S+QnVzwSFeW5XLgtW5/M/CyBdaJsQZI9I6cV7/rpyX3pWRaZ1ienmrvDxyBjRn2Q7eXJfHsTLn0sHd+MoFp5PRu7MGvqVB0ZmHSC32HS7h4x37yNy+j79v/ozV2fspd2jfKoELBiQFEzR2o1Ob8DP9lpc7K3buZ/6qXSxYnUvegWI6t2nBtWemMnVMGv27tT8JPRKJns48RE5Q57aJXDK4+z+m8SgsOsaHW/bw1oYC/rYxn/mrckmIM85L78qVw3sytm8Xktq3/MedW2Xlzv6iEnL2H2FzwSG2FByOvPYcZtuewxw5VkZifBwXDEziiyN6ctnQ7hoAlwZPZx4idVBe7qzOKWTB6lzmr8olZ/+Rf2zr1KYFBuw/coyK/83i44y0zq3pl9SOvl3bMiylIxcP7kaHVhoAl4ZHZx4iJ0FcMA4yIq0TMyYO4pPsQj7NO0j+gaPsPnAUdzitbSJd2ibSo2Nr+ndrS68ubZvlHFvStCg8RGLEzBiZ1omRaZ3quykiJ51+/RERkdAUHiIiEprCQ0REQqtTeJhZFzN708w2Be+dqyk3wcw2mlmWmc2orb6Z3WJmKyu8ys1sZLDt7WBfx7c1/S9nEBFpYOp65jEDWOzu6cDiYPlfmFk88AgwERgCTDWzITXVd/dn3X2ku48EbgO2ufvKCru95fh2d8+vYx9ERCSkuobHZGBW8HkWcHUVZcYAWe6+xd1LgDlBvWjrTwWer2M7RUQkhuoaHt3dPRcgeK/qElIKsLPCcnawLtr6N/H58Ph9cMnqAdOEPyIip1ytz3mY2SKgRxWbvh/lMar64R7VY+1mNhYocvc1FVbf4u45ZtYe+BORy1qzq6k/HZgO0KtXryibKyIitak1PNz90uq2mVmemSW7e66ZJQNVjT9kA2kVllOBXcHn2upPodJZh7vnBO8Hzew5IpfFqgwPd58JzITI9CTV9UNERMKp62WrecAdwec7gLlVlFkGpJtZXzNLJBII82qrb2ZxwA1ExkiOr0sws67B5xbAlUDFsxIRETkF6hoeDwLjzWwTMD5Yxsx6mtkCAHcvBe4FFgLrgRfdfW1N9QPjgGx331JhXUtgoZmtAlYCOcATdeyDiIiEpFl1RUSkWtXNqqsnzEVEJDSFh4iIhKbwEBGR0BQeIiISmsJDRERCU3iIiEhoCg8REQlN4SEiIqEpPEREJDSFh4iIhKbwEBGR0BQeIiISmsJDRERCU3iIiEhoCg8REQlN4SEiIqEpPEREJDSFh4iIhKbwEBGR0OoUHmbWxczeNLNNwXvnaspNMLONZpZlZjMqrL/BzNaaWbmZZVSqc39QfqOZXV5h/WgzWx1se9jMrC59EBGR8Op65jEDWOzu6cDiYPlfmFk88AgwERgCTDWzIcHmNcC1wLuV6gwBpgBDgQnA74L9ADwKTAfSg9eEOvZBRERCqmt4TAZmBZ9nAVdXUWYMkOXuW9y9BJgT1MPd17v7xmr2O8fdi919K5AFjDGzZKCDu3/o7g7MruaYIiJyEiXUsX53d88FcPdcM+tWRZkUYGeF5WxgbC37TQGWVKqTAhwLPldeXyUzm07kLAXgkJlVFVTR6ArsOcG6jZX63Dw0tz43t/5C3fvcu6qVtYaHmS0CelSx6ftRHriqMQk/wTqh9uXuM4GZtRyrVmaW6e4ZtZdsOtTn5qG59bm59RdOXp9rDQ93v7S6bWaWZ2bJwVlHMpBfRbFsIK3Cciqwq5bDVlcnO/gcZl8iIhJjdR3zmAfcEXy+A5hbRZllQLqZ9TWzRCID4fOi2O8UM2tpZn2JDIwvDS6RHTSzs4O7rG6v5pgiInIS1TU8HgTGm9kmYHywjJn1NLMFAO5eCtwLLATWAy+6+9qg3DVmlg2cA7xmZguDOmuBF4F1wF+Be9y9LDjmV4EniQyibwZer2MfolHnS1+NkPrcPDS3Pje3/sJJ6rNFbloSERGJnp4wFxGR0BQeIiISmsKjBtVNq9KUmFmamb1lZuuDqWLuC9ZHNfVMY2Zm8Wa2wszmB8tNus9m1snMXjazDcHf9znNoM/fCv5drzGz582sVVPrs5k9bWb5Zramwrpq+1jd1E9hKTyqUcu0Kk1JKfB/3H0wcDZwT9DPWqeeaQLuI3ITx3FNvc8PAX9190HACCJ9b7J9NrMU4BtAhrufAcQTuduzqfX5GT4/TVOVfaxl6qdQFB7Vq3ZalabE3XPd/ePg80EiP1BSiG7qmUbLzFKBK4jcuXdck+2zmXUAxgFPAbh7ibvvpwn3OZAAtDazBKANkefCmlSf3f1dYG+l1dX1scqpn07kuAqP6lU1rUq1U6E0BWbWBxgFfESlqWeAqqaeacx+DXwXKK+wrin3uR9QAPw+uFT3pJm1pQn32d1zgF8AO4BcoNDd36AJ97mC6voYs59rCo/qnci0Ko2WmbUD/gR8090P1Hd7TiYzuxLId/fl9d2WUygBOBN41N1HAYdp/JdrahRc558M9AV6Am3N7Nb6bVW9i9nPNYVH9U5kWpVGycxaEAmOZ939lWB1XjDlDDVMPdNYnQt80cy2EbkcebGZ/ZGm3edsINvdPwqWXyYSJk25z5cCW929wN2PAa8AX6Bp9/m46voYs59rCo/qnci0Ko1OMM3LU8B6d/9VhU3RTD3TKLn7/e6e6u59iPy9/s3db6Vp93k3sNPMBgarLiEyg0OT7TORy1Vnm1mb4N/5JUTG9Jpyn4+rro9VTv10IgfQE+Y1MLNJRK6NxwNPu/tP67dFsWdm5wHvAav55/X//yAy7vEi0IvIf8Ib3L3yoFyjZ2YXAt9x9yvN7DSacJ/NbCSRGwQSgS3ANCK/QDblPv8ncBORuwpXAHcB7WhCfTaz54ELiUy9ngf8EPgL1fTRzL4PfInIn8k33f2EpnhSeIiISGi6bCUiIqEpPEREJDSFh4iIhKbwEBGR0BQeIiISmsJDRERCU3iIiEho/x8sBhMqSrl4iQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "def animate(i):\n", " ax.clear()\n", " ax.plot(sol[i*10])\n", " ax.set_ylim(-0.01, 0.01)\n", " \n", "fig, ax = plt.subplots(1,1)\n", "ax.set_ylim(-0.01, 0.01)\n", "ani = animation.FuncAnimation(fig, animate, frames=500, interval=50)\n", "ani.save('string.gif',writer='pillow',fps=20)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Create .WAV file of noise" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Extract the \"amount\" of the harmonics at any time $t$:\n", "\n", "$$ \\text{Amplitude of harmonic n at time t} \\propto \\int_{0}^L y(x, t) \\sin(n \\pi x / L) dx $$" ] }, { "cell_type": "code", "execution_count": 119, "metadata": {}, "outputs": [], "source": [ "# This is the way I did it in the video\n", "def get_integral(n):\n", " sin_arr = np.sin(n*np.pi*np.linspace(0,1,101))\n", " return np.array([sum(sin_arr*s) for s in sol])\n", "\n", "# This is the same as the function above. but runs WAAYYY faster. Use this instead; it is far more optimized\n", "def get_integral_fast(n):\n", " sin_arr = np.sin(n*np.pi*np.linspace(0,1,101))\n", " return np.multiply(sol, sin_arr).sum(axis=1)" ] }, { "cell_type": "code", "execution_count": 120, "metadata": {}, "outputs": [], "source": [ "hms = [get_integral_fast(n) for n in range(10)]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Add them together" ] }, { "cell_type": "code", "execution_count": 126, "metadata": {}, "outputs": [], "source": [ "all_harmonics=True\n", "if all_harmonics:\n", " tot = sol.sum(axis=1)[::10] # all harmonics\n", "else:\n", " tot = sum(hms)[::10] # only first 10 harmonics\n", "tot = tot.astype(np.float32)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Make a WAV file" ] }, { "cell_type": "code", "execution_count": 127, "metadata": {}, "outputs": [], "source": [ "from scipy.io import wavfile\n", "from IPython.display import Audio" ] }, { "cell_type": "code", "execution_count": 128, "metadata": {}, "outputs": [], "source": [ "wavfile.write('sound.wav',20000,tot)" ] }, { "cell_type": "code", "execution_count": 129, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "execution_count": 129, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Audio('sound.wav')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Only on my computer: You will have to create these files" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "from scipy.io.wavfile import read\n", "a = read('../output/v3_Al.wav')[1]\n", "c = read('../output/v3_Csl.wav')[1]\n", "f = read('../output/v3_Fsl.wav')[1]\n", "d = read('../output/v3_D.wav')[1]" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "wavfile.write('../output/v3_Al+Csl.wav',20000,a+2*c)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Audio('../output/v3_Al+Csl.wav')" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "a = np.append(a, np.zeros(14000))\n", "c = np.concatenate([np.zeros(1000), c, np.zeros(14000-1000)])\n", "d = np.append(np.zeros(14000), d)\n", "f = np.concatenate([np.zeros(14000-1000), f, np.zeros(1000)])\n", "tot = a + c + d + 0.4*f\n", "tot = tot.astype(np.float32)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "wavfile.write('../output/v3_Al+Csl+D.wav',20000,tot)" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Audio('../output/v3_Al+Csl+D.wav')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Alternative Solution" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "IN PROGRESS: NOT WORKING YET" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As recommended by https://www.reddit.com/r/Physics/comments/miabsi/solving_the_full_damping_stressstrain_wave/gt434hj?utm_source=share&utm_medium=web2x&context=3" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy.integrate import odeint" ] }, { "cell_type": "code", "execution_count": 107, "metadata": {}, "outputs": [], "source": [ "L = 0.7\n", "f = 220\n", "c = 2*L*f \n", "l=2e-6\n", "gamma=2.6e-5\n", "N = 101\n", "dx = L/(N+1)\n", "x = np.arange(0, L, N)" ] }, { "cell_type": "code", "execution_count": 108, "metadata": {}, "outputs": [], "source": [ "l = 0\n", "gamma = 0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Wave numbers" ] }, { "cell_type": "code", "execution_count": 109, "metadata": {}, "outputs": [], "source": [ "kappa = 2 * np.pi * np.fft.fftfreq(N, dx)" ] }, { "cell_type": "code", "execution_count": 110, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "4.487989505128276" ] }, "execution_count": 110, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.pi / L" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Initial Condition" ] }, { "cell_type": "code", "execution_count": 139, "metadata": {}, "outputs": [], "source": [ "ya = np.linspace(0, 0.01, 30)\n", "yb = np.linspace(0.01, 0, 21)\n", "y0 = np.concatenate([ya, yb[:-1], np.array([0]), -np.flip(yb)[1:], -np.flip(ya)])\n", "#y0 = 1e-3 * np.sin(np.pi*np.linspace(-1, 1, N))\n", "y0_h = np.fft.fft(y0)" ] }, { "cell_type": "code", "execution_count": 140, "metadata": {}, "outputs": [], "source": [ "y0_h_ri = np.concatenate((y0_h.real, y0_h.imag))\n", "z0_h_ri = np.zeros(len(y0_h_ri))\n", "S0 = np.concatenate([y0_h_ri, z0_h_ri])" ] }, { "cell_type": "code", "execution_count": 143, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "C:\\Users\\lukep\\anaconda3\\lib\\site-packages\\numpy\\core\\_asarray.py:83: ComplexWarning: Casting complex values to real discards the imaginary part\n", " return array(a, dtype, copy=False, order=order)\n" ] }, { "data": { "text/plain": [ "[]" ] }, "execution_count": 143, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAD4CAYAAAAUymoqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA5gUlEQVR4nO3dd3gV1dbA4d9KpfcQelMsoBDCSQAV9CooCEpvgjQVUVEpesWreO1dvNcCCBZAmlgoIlLE3gmK0iX0UEIAU4CQur4/Mnw3xoQk5CSTst7nOc85M7P3zNqUrJwpa4uqYowxxuSHj9sBGGOMKXkseRhjjMk3Sx7GGGPyzZKHMcaYfLPkYYwxJt/83A6gqNSqVUubNGnidhjGGFOirF+//qiqBmVdX2aSR5MmTYiIiHA7DGOMKVFEZG926+20lTHGmHyz5GGMMSbfLHkYY4zJN0sexhhj8s2ShzHGmHzzSvIQka4isl1EIkVkUjbbRURecbb/LiKhmba9LSJHRGRTlj41RGSNiOxw3qtn2vags6/tInKdN8ZgjDEm7wqcPETEF3gd6Aa0AAaLSIsszboBzZ3XaGBapm2zgK7Z7HoSsFZVmwNrnWWcfQ8CWjr9pjoxGGOMKSLeeM4jHIhU1V0AIrIQ6AlsydSmJzBHM+q//ygi1USkrqoeUtWvRaRJNvvtCVzlfJ4NfAk84KxfqKpJwG4RiXRi+MELYzFFLWY7bPoQVDkcf5pdPo1p3+MWfHzE7ciMMWfhjeRRH9ifaTkKaJeHNvWBQ2fZb7CqHgJQ1UMiUjvTvn7MZl9/IyKjyfimQ6NGjc4+ClP0ju2Ed66HU0dJR6iDUgeYcjCV0SNvpVJgmXmG1ZgSxxvXPLL7FTHrDFN5aePN42WsVJ2hqh5V9QQF/e3peuOmE0dgbh9Op6RxbeoUWqQt4L/tvyGufEP6HHqZga9/wd5jJ92O0hiTA28kjyigYablBsDBc2iTVbSI1AVw3o8UYF+mOEk6AfMHkBofzeCT46nT7FK+uO8q7u3aiqr9XqGJHKZHwiJ6T/2emIQkt6M1xmTDG8ljHdBcRJqKSAAZF7OXZWmzDBjm3HXVHog7c0rqLJYBw53Pw4GlmdYPEpFAEWlKxkX4n70wDlMU0lJg0TD00O/cmXQ3vg3DmHFzW+pWLZ+x/byr4ZK+3C5LqHF6H0+v2OpuvMaYbBU4eahqKjAWWAVsBRap6mYRGSMiY5xmK4BdQCQwE7jzTH8RWUDGxe4LRSRKRG5xNj0LdBGRHUAXZxlV3QwsIuOC/ErgLlVNK+g4TBFQhWV3w861PKq3srP65bw53EM5/yw3y133ND7+5Xir1kIW/xrF9zuPuhOvMSZHknEDVOnn8XjUquq67LPH4NspzPAdxFu+/fnwjstoUL1C9m1/ngkr7uPRgIl8HdiJT+/tSKCf3ZFtTFETkfWq6sm63p4wN0Xj55nw7RRW+F/Lq6m9eWdEeM6JA8AzCuq14UGfOcTExDDjq11FF6sxJleWPEzh27IMXXE/EYHtmHBqGG/c7KFFvSpn7+PjCz1eJjDpGK8GL+fVLyLt7itjihFLHqZw7f0B/fBWdpe7mKFxY3h+QFsuO79W3vrWawNht3Fl3FJCfHbxyNLNlJXTrMYUd5Y8TOE5sg1dMIhjfsH0jb2H+7qHcGPrevnbx9UPIZWCmVr1Xb75I5oVGw8XTqzGmHyx5GEKR/xBmNuXU+l+9IqfQN8rWnNrx2b530+5qtD1aWolbOW+6t/w+PLNJJxO8X68xph8seRhvC8xFub2I+XkcQYkTKBNqxD+df3F576/ln3gvKu5PXUeJBxmypo/vBaqMebcWPIw3pWaBO8NJT1mO7eevpcqTdvyYv9WBSt0KALXv4hvegozgz9i9vd72HQgznsxG2PyzZKH8Z70dFh8O+z5hklpY4gOuow3hrX1zvMZNc+DjhNpFbuWbuW38NDijaSl28VzY9xiycN4hyqsfgg2L+a/PjfzXYVrmD0qnCrl/L13jCvGQc3zea78bLZFxTD/533e27cxJl8seRjv+OE1+HEqH/jfwNvpNzB7VBjBVcp59xh+gdD9JSqd3M/TQWt4fuU2jiSc9u4xjDF5YsnDFNzGD2D1w3wbeAUPJw7mrRFhnF+7cuEcq9lVcGl/+px8n7opUTz1iRVONMYNljxMwez6El08hm2Brbg1/lb+O7gtniY1CveY1z6F+JdnZq0FLN1wgO8irXCiMUXNkoc5d4d+RxcOJdq/IQPi7ubhnqFc17JO4R+3cjBcM5nGceu4pWoEk5dsIinVCisbU5QseZhz8+demNePE5SnV9x4hv2jNUPbNy6643tGQb1QHpB3OXr0CG9Y4URjipQlD5N/p47D3L4knz5Fn4T76Ni2NROvvaBoY3AKJwYkHefV4OW89kUke45a4URjioolD5M/yadg/kDS/tzLzafGU/+CNjzd51JECvAQ4LmqFwLho+kUt4xQ311MXrrJCicaU0S8kjxEpKuIbBeRSBGZlM12EZFXnO2/i0hobn1F5D0R2eC89ojIBmd9ExFJzLRtujfGYPIgLRU+vAWNWse4lLGcrteOqUNC8fd18XeQf2QUTny9yhy+3xHNJxtzm93YGOMNfgXdgYj4Aq+TMVVsFLBORJap6pZMzbqRMdd4c6AdMA1od7a+qjow0zFeAjLXo9ipqiEFjd3kgyqsuA+2r+A5GcXGKp34YEQYFQIK/E+oYMpVgW7PUvP9ETxQ41se/7gCV14QRGVvPpxojPkbb/zKGA5EquouVU0GFgI9s7TpCczRDD8C1USkbl76Ssb5kAHAAi/Eas7V1y/C+nd4168P7/tcz6yR4dSqFOh2VBla9ILzO3NLyjx8ThzipdVWONGYwuaN5FEf2J9pOcpZl5c2eenbEYhW1R2Z1jUVkV9F5CsR6ZhTYCIyWkQiRCQiJiYmb6Mxf/fLu/DFk3wWcDXPJA/gnZFhNKlV0e2o/kcErn8BX03lzeCPmPPDHjZGWeFEYwqTN5JHdldKs161zKlNXvoO5q/fOg4BjVS1DTABmC8i2c5pqqozVNWjqp6goKBsgze5+GM1+vG9/B7YlrtOjGTqkLa0alDN7aj+rkYz6Hgfl8R+To8KW3hoiRVONKYweSN5RAENMy03AA7msc1Z+4qIH9AHeO/MOlVNUtVjzuf1wE6giO8TLSOi1qPvD2d/4HkMjruTp/qGctWFtd2OKmeX3wM1m/NM4Gy2R8Uw76e9bkdkTKnljeSxDmguIk1FJAAYBCzL0mYZMMy566o9EKeqh/LQtzOwTVWjzqwQkSDnQjsi0oyMi/D2hJi3HdsJ8/sT61OdPrHjufO6EPq1beB2VGfnFE6seGo/zwWt5oWV261wojGFpMDJQ1VTgbHAKmArsEhVN4vIGBEZ4zRbQcYP+EhgJnDn2fpm2v0g/n6hvBPwu4j8BnwAjFHV4wUdh8nkxBGY24fElHR6x0+gW/tW3HnVeW5HlTfNroRWA+l58n3qp+3nyeVWONGYwiBl5aEqj8ejERERbodR/CUlwKzupB75g36JDxJ88eVMHdIW34LMBFjUThyB1zzsC2xOp+jxvHtLOzo2t2texpwLEVmvqp6s6+0Jc/M/aSmwaDh6eBN3Jt2NX8Mw/juoTclKHACVasM1/6ZRXAS3VV3HI0s3czrFCica402WPEwGVVh2N+xcyyN6G7tqXMGbwz2U8/fCFLJuaDsS6nu4nzkcOxrN9K92uh2RMaWKJQ+TYe1j8NsC3vAdxOqALswaGUa1CgFuR3XufHwyCicmx/J68HKmfrmT3VY40RivseRh4KcZ8O3LLPe/jtdSezNrZDgNqldwO6qCq9sK2o3hiriPCfON5BErnGiM11jyKOu2LEU//SfrAttz/6lhvDHMw8V1s33msmT6x7+QynV5tcq7fL8jmo9/t8KJxniDJY+ybO/36Ie3savcxdwcdzvPDwjlsvNquR2VdwVWhq7PUCNhOw/W/Jonlm8h/nSK21EZU+JZ8iirjmxFFwziqF8wfWPv5f4ebbihdT23oyocLXrC+V0YmTwf/xMHeWnVdrcjMqbEs+RRFsUdgLl9OZXuR+/4ifTv2IpbrmjqdlSF5/8LJ6YxM/gj5vy4l9+jYt2OypgSzZJHWZMYC/P6kXIqlv4JEwlt1ZoHu13sdlSFr0ZT6HQfLWO/4MYKm3lo8SYrnGhMAVjyKEtSTsPCIaQf3cGoxHFUP68tL/ZvjU9JewjwXF12L9S6gKcDZ/PHgRjm/miFE405V5Y8yor0dFh8O+z9lgdSx3C0dgemD21LgF8Z+ifgFwDdp1DxVBTP117FC6u2Ex1vhRONORdl6CdHGaYKq/4FW5bwss9wvq9wNbNHhpXNqVqbdoRWg7jxxAc0SNvPk59Y4URjzoUlj7Lg+1fhp2m8738Ds7U7s0eFU7tKObejcs+1TyIBFZhZcz4f/3aAr/+wWSaNyS9LHqXd74tgzWS+DezI5MTBvDUijPNrV3I7KndVCoLOj9Ewbj2jq67jkaWbrHCiMflkyaM02/kFuuROtpZrza3xt/LK4La0bVzD7aiKh9Dh0CCM+5jD8WNHmPalFU40Jj8seZRWh35D3xvKYf+GDIwdy8M923BtyzpuR1V8/H/hxDimBn/MtC93sivmhNtRGVNieCV5iEhXEdkuIpEiMimb7SIirzjbfxeR0Nz6isijInJARDY4r+szbXvQab9dRK7zxhhKlT/3wLz+JFCRXnETGHF1a4a2b+x2VMVPnUuh/R0ZhRP9I5lshRONybMCJw9nPvHXgW5AC2CwiLTI0qwbGXONNwdGA9Py2PdlVQ1xXiucPi3ImJ62JdAVmHpmTnMDnDwGc/uSlJRIn4SJXOlpxfguF7gdVfF11SSoXI9XK7/Lj5FHWPbbQbcjMqZE8MY3j3AgUlV3qWoysBDomaVNT2COZvgRqCYidfPYN6uewEJVTVLV3WTMix7uhXGUfMmnYMFA0v7cx80nx9HwgjY81ftSRMrIQ4DnIrAydHuOGgnb+VfNr3hi+VbiEq1wojG58UbyqA/sz7Qc5azLS5vc+o51TnO9LSLV83E8AERktIhEiEhETEwpvx0zLRU+vAWNimBcyliS6rXj9SGh+PvaZa1cXXwDNL+OEUnzCTx5kBetcKIxufLGT5bsfq3NeuI4pzZn6zsNOA8IAQ4BL+XjeBkrVWeoqkdVPUFBQdk1KR1U4ZMJsH0Fz8ooNlbpxNsjwqgQ4Od2ZCWDCFz/PL4oM4M/ZO5Pe9mwP9btqIwp1ryRPKKAhpmWGwBZTxzn1CbHvqoarappqpoOzOR/p6bycryy5avn4ZfZzPHty4e+3Zg9KpyalQLdjqpkqd4ErryfFrFf0rviJh5avJHUtHS3ozKm2PJG8lgHNBeRpiISQMbF7GVZ2iwDhjl3XbUH4lT10Nn6OtdEzugNbMq0r0EiEigiTcm4CP+zF8ZRMq2fDV8+zZqAa3g2pT9vjwijcc2KbkdVMnW4G2pdyJMBs9l5MIZ3rXCiMTkq8HkNVU0VkbHAKsAXeFtVN4vIGGf7dGAFcD0ZF7dPASPP1tfZ9fMiEkLGKak9wO1On80isgjYAqQCd6lq2Xw8+I9V6PLx/Bbo4e6EEbwxwkOrBtXcjqrk8guAHi9TYdb1vFB7FQ+urki3S+pSp2oZLuViTA6krNzX7vF4NCIiwu0wvCcqAp3Vg/2+DekaN4kn+rejb9sGbkdVOiy5E/39PbonP0PTiz28PiQ09z7GlFIisl5VPVnX2604JdHRSJg/gD99a9Inbjxju7a2xOFNXR5HAioxs8YCPtl4kC+3H3E7ImOKHUseJU1CNMztQ2JKOr3jJ9C9QyvuuPI8t6MqXSrWgi6PUz/+F8ZU+5lHlm62wonGZGHJoyRJSoD5/UlNOMKgE+O5uEUIj9zQ0h4CLAxtboaG7Zio7xJ/PJqpX0S6HZExxYolj5IiNRkWDUMPb2JM0j0ENA7jP4NC8C0rU8gWNadwon9yHNPqLGPaVzuJPGKFE405w5JHSaAKy+6GnZ/ziI5mb43LeXNYGOX8raRXoQpuCR3upEPsJ7T3j2TyEiucaMwZljxKgrWPwe8Lme47mDUBXZg9KpyqFcrgFLJuuHISVGnAK5XmsG5XNEs2HHA7ImOKBUsexd1PM+Dbl/nYvyuvp/Vi1qgw6lUr73ZUZUdgJej2HNVP7OBfNb/mqU+2EnfKCicaY8mjONuyFP30n/wc2IH7Tw1jxs1hXFSnittRlT0XdYcLujE8aT7lTh7khdXb3I7IGNdZ8iiu9n6Pfngbu8q1YFj87bw4sA0dzqvpdlRl05nCiQJvBn/AvJ/28eu+P92OyhhXWfIojo5sRRcM4qhfMH1j7+Gf3UPo0aqe21GVbdUawZX/5KLYr+lT8XceWrzJCieaMs2SR3ETFwVz+3Iq3Z9ecRMZ0Kk1o65o6nZUBqDDWAi6mCcC5rD70BFm/2CFE03ZZcmjOEn8E+b2I+VULP0SJhIW0ppJXS9yOypzhq8/9JhChVMHebH2Kqas3s7huNNuR2WMKyx5FBcpp2HhENKPRTIqcRw1zgvl+X6t8bGHAIuXxpdByFCuP/EhTdL38fjyzbn3MaYUsuRRHKSnweLRsPc77k+9g2O1OzB9aFsC/Oyvp1jq8jgSWJmZNebx6caDfGGFE00ZZD+d3KYKKx+ELUuZIsP5scI/mDUyjMrl7CHAYqtiTejyOPXiN3BHtZ/4txVONGWQV5KHiHQVke0iEikik7LZLiLyirP9dxEJza2viLwgItuc9otFpJqzvomIJIrIBuc13RtjcM13/4Wf32CR343MkR7MuSWc2lVs8qFiL2QoNGzPeJ1L/PFoXvvcCieasqXAyUNEfIHXgW5AC2CwiLTI0qwbGdPFNgdGA9Py0HcNcImqtgL+AB7MtL+dqhrivMYUdAyu+W0hfPZvvg7oxL+TBvHW8DDOC6rkdlQmL3x8oMcU/JPjmR68jDe+tsKJpmzxxjePcCBSVXepajKwEOiZpU1PYI5m+BGo5sxRnmNfVV2tqqlO/x+B0jXb0c7P0aV3sSUwhNEJt/DKYA9tG1d3OyqTH8EtocNdtI/7hMv8d/Dwko1WONGUGd5IHvWB/ZmWo5x1eWmTl74Ao4BPMy03FZFfReQrEemYU2AiMlpEIkQkIiYmJveRFJWDG9D3buaQf2MGxo3lkV6hdGkR7HZU5lxclVE48b+V5hCx64gVTjRlhjeSR3b3kmb99SunNrn2FZGHgFRgnrPqENBIVdsAE4D5IpJtwSdVnaGqHlX1BAUFnWUIRejPPTCvPwlUolfcBEZe05qb2jVyOypzrgIqwvXPU+1EJJNrfsmTy61woikbvJE8ooCGmZYbAAfz2OasfUVkONADGKLO+QBVTVLVY87n9cBO4AIvjKPwnTwG7/YhKfk0vRMmcnVYK8Z3bu52VKagLuoOF17PzUkLKJ94iOdWWeFEU/p5I3msA5qLSFMRCQAGAcuytFkGDHPuumoPxKnqobP1FZGuwAPAjap66syORCTIudCOiDQj4yL8Li+Mo3Aln4L5A0iLjWLoyfE0vrANT/a6xKaQLS26PYePCG/Vfp8FP+/jFyucaEq5AicP56L2WGAVsBVYpKqbRWSMiJy5E2oFGT/gI4GZwJ1n6+v0eQ2oDKzJcktuJ+B3EfkN+AAYo6rHCzqOQpWWCh+MQg/+wt0pY0mp347XbmqDn689ZlNqVGsEVz7AhbFf088KJ5oyQMrK3SEej0cjIiKK/sCq8PE98MscnuRW1la6gQ/GdKBmpcCij8UUrrQUmN6RxJNxhB5/kond23Brx2ZuR2VMgYjIelX1ZF1vv/oWtq+eg1/mMMu3H0v8ujF7ZLgljtLK1x96vEz5Uwd5qfYqpqz5g4OxiW5HZUyhsORRmNbPhi+fYXXANbyQ0o9ZI8NoVLOC21GZwtS4A7QZSrcTH3K+7uWxj61woimdLHkUlu0r0eXj2BDo4Z4TI5g21MMl9au6HZUpCp0fRwKrML36PFZvPsTn26LdjsgYr7PkURj2r0PfH8G+wObcFHcnT/cLpdMFxeQ5E1P4KtaEa5+gXvxvjK32I48s3UxishVONKWLJQ9vOxoJ8wfwp28N+sSO4+6uIfQJLV2VVUwetL4JGl3GPelzOflnNK9+vsPtiIzxKkse3pQQDXN7k5iq9I6fyA2XtWbMlXa3TZl0pnBi6gneqLOMGV/vYkd0gttRGeM1ljy85XQ8zOtHasJRBpyYSIuWIUzu0cIeAizLal8MHcYSHruCKwJ28NCSTVY40ZQaljy8ITUZFt2MRm9mTNLdlG/s4eWBIfjaFLLmyn9C1Ub8p9Icft19hA9/scKJpnSw5FFQ6emw9C7Y9SWT029nX83LmTnMQzl/X7cjM8VBQEW4/gWqnYjkkZpf8PSKrcSeSnY7KmMKzJJHQa19FDYuYprPYD4L7MyskeFUrWBTyJpMLuwKF/VgSNICKiUe5LmVVjjRlHyWPArix+nw3X9Z6t+Nqem9mD0qnHrVyrsdlSmOuj6Lj/jyVtB7LPh5H+v3Fu9ybMbkxpLHudq8BF05iZ8CO/DAqZuZOSyMC+tUdjsqU1xVawhXTaJ53HcMrJRRODHFCieaEsySx7nY8y360W3sDGzJ8PjbeWlgKO2b1XQ7KlPctb8DarfkUf857Dscw6zv9rgdkTHnzJJHfkVvQRcMJsavHn3j7uGBHiF0b1XX7ahMSXCmcGLiIV6u/Skvf2aFE03JZckjP+KiYG5fTmogveImMKhTK0Ze3tTtqExJ0qgdhA7j2oSPOF/3WOFEU2JZ8sirxD9hbl+SE+PplzCR8JBWPND1IrejMiVR58eQ8tV4o/p8Vm8+xGdbrHCiKXm8kjxEpKuIbBeRSBGZlM12EZFXnO2/i0hobn1FpIaIrBGRHc579UzbHnTabxeR67wxhrNKOQ0LbiL92E5GJo4j6PxQnu/XGh97CNCciwo14NonqRv/G3dX+5F/L9vMqeRUt6MyJl8KnDyc+cRfB7oBLYDBItIiS7NuZMw13hwYDUzLQ99JwFpVbQ6sdZZxtg8CWgJdgaln5jQvFOlpnFw4CvZ9zz2nxxBXpz3ThrYlwM++tJkCaD0YGl/OPTqXxNhoXv080u2IjMkXb/wEDAciVXWXqiYDC4GeWdr0BOZohh+BaiJSN5e+PYHZzufZQK9M6xeqapKq7iZjXvRwL4zjb1JT01g3fTQVd37Cs+nDaH71cBbd3oFKgX6FcThTlohA9yn4pZzgjTpLmPn1Lv6wwonG2/Z8BwuHQNIJr+/aG8mjPrA/03KUsy4vbc7WN1hVDwE477XzcTwARGS0iESISERMTEyeB3SGH2kEJOznq5oDGXnfC9zbuTkVAixxGC+pfRFcdg9hsSvpFPgHDy+2wonGi6K3oAsGkXZkO6R5vySON5JHdif+s/4PyKlNXvqey/EyVqrOUFWPqnqCgs5hMia/AFpOWM6Vd00nuEq5/Pc3Jjed7odqjXi54hx+3XOED9ZHuR2RKQ3iomBeP05pAL3iJ3Aw2fuVL7yRPKKAhpmWGwAH89jmbH2jnVNbOO9H8nE8r/HzD8iYm8GYwhBQAa5/kaondvLvWhmFE/88aYUTTQEkxsLcfqSciqNfwkTOO/9i6hTCL7/e+Km4DmguIk1FJICMi9nLsrRZBgxz7rpqD8Q5p6LO1ncZMNz5PBxYmmn9IBEJFJGmZFyE/9kL4zDGHRdcBxffwE2nF1Dl9EGe/dQKJ5pzlHIaFt5E+rFIRibeS83z2hbanaEFTh6qmgqMBVYBW4FFqrpZRMaIyBin2QpgFxkXt2cCd56tr9PnWaCLiOwAujjLONsXAVuAlcBdqmoTRJuSzSmc+GbQIt6L2EfEHiucaPIpPR0Wj4a933F/6h0cr92BaUNDC+3OUCkrF+g8Ho9GRES4HYYxOfv+NVj9EJP8/smvFTuy/J4r8Pe1U6YmD1Rh5ST4aTov+wzng4BeLL7zMmp74XSViKxXVU/W9fYv05jiot0YCL6UR/3nsD86hne+2+12RKak+P4V+Gk67/ndyGx6MOeWcK8kjrOx5GFMceHrBz2mUC7xcEbhxDU7OGCFE01ufnsP1jzCNwGdeDRpEG8ND+O8oEqFflhLHsYUJw3Doe0Irk34iAvZw6PLrHCiOYudn6NL72RrYAi3JdzCK4M9tG1cPfd+XmDJw5ji5pp/I+WrM736XD7bcog1VjjRZOfgBvS9mzkU0JgBcWN5pFcoXVoEF9nhLXkYU9w4hRPrxG/knuo/8KgVTjRZHd8N8/qTQEV6xU5g5DWtualdoyINwZKHMcVR60HQ+AruTpvL6djD/HftDrcjMsXFyWMwty9JyafpnXAf//C0Ynzn5kUehiUPY4ojEegxBb/UU8yos5Q3v9nN9sNWOLHMSz4F8weQFhvF0JPjaXxhG57qfQkiRT89hCUPY4qroAvh8ntpG7uSqwO38fCSjaSnl43nskw20lLhg5HowV+4O2UsyfXb8dpNbfBz6VkgSx7GFGed7oPqTXip4hw27ImxwolllSp8Mh7+WMlTOoqtVTvx9nCPq1W+LXkYU5z5l4frX6TKid08Xutznv50K8etcGLZ89Vz8MscZvn2Y4lfV2aPDKdmpUBXQ7LkYUxx17wLtOjJwMQFVD99gGdWbHU7IlOU1s+CL59hdcA1vJDSj3dGhNOoZgW3o7LkYUyJ0PVZfHz9eTPoPd5fv5+fd1vhxDJh+0p0+Xg2BHq458QIpg31cGmDqm5HBVjyMKZkqFIP/vEQ58X9wJDKv/Hwko2kpKW7HZUpTPvXoe+PYF9gc26Ku5On+4XS6YJzmNSukFjyMKakCB8NdS5lst8sDkYf4a1vrXBiqXU0EuYP4E/fmvSJHcfdXUPoE9rA7aj+wpKHMSWFrx/0+A/lEmP4T+1P+c9nf7D/+Cm3ozLelhANc3uTmKr0ip/IDZe1ZsyVzdyO6m8seRhTkjTwgGck1yQspgV7eOxjK5xYqpyOh3n9SE04yoATE7nkktZM7tHClYcAc1Og5CEiNURkjYjscN6zLecoIl1FZLuIRIrIpNz6i0gXEVkvIhud96sz9fnS2dcG51W7IGMwpsS55hGkQk2mV5/L51sPs3rzYbcjMt6QmgyLbkajN3N70j2Ub+xhyoAQfAthCllvKOg3j0nAWlVtDqx1lv9CRHyB14FuQAtgsIi0yKX/UeAGVb2UjPnL382y2yGqGuK8jhRwDMaULOWrw7VPUTt+E+Oqf8+jyzZzMskKJ5Zo6emwbCzs+pLJ6bezv+ZlzBzmoZy/r9uR5aigyaMnMNv5PBvolU2bcCBSVXepajKw0OmXY39V/VVVDzrrNwPlRMTdJ2KMKU5aDYAmHbkzbS7JcdG8YoUTS7a1j8Lv7zHNZzCfBXZm1shwqlbwdzuqsypo8ghW1UMAznt2p5DqA/szLUc56/Lavy/wq6omZVr3jnPKarKc5WSgiIwWkQgRiYiJicn7qIwp7kSg+xT8UhOZUWcxb367m22H492OypyLH6fDd/9lqf/1TE3vxexR4dSrVt7tqHKVa/IQkc9EZFM2r5659T2zi2zW5am6m4i0BJ4Dbs+0eohzOquj87o5p/6qOkNVParqCQoqPvdHG+MVQRfAFeMIjV1N53LbeGjxJiucWNJsXoyunMRPgZfxwKmhzBwWxoV1KrsdVZ7kmjxUtbOqXpLNaykQLSJ1AZz37K4/RAENMy03AM6cksqxv4g0ABYDw1R1Z6Z4DjjvCcB8Mk6LGVM2dZwI1ZvwQvk5bNx7hPfX78+9jyke9nyLfjSanYEtGR4/mpcGhtK+WU23o8qzgp62WkbGBW2c96XZtFkHNBeRpiISAAxy+uXYX0SqAZ8AD6rqd2d2JCJ+IlLL+ewP9AA2FXAMxpRc/uXh+peocnI3j9dayzOfbuPYiaTc+xl3RW9GFwwmxq8efePu4YEeIXRvVdftqPKloMnjWaCLiOwAujjLiEg9EVkBoKqpwFhgFbAVWKSqm8/W32l/PjA5yy25gcAqEfkd2AAcAGYWcAzGlGzNO0OLXgxIfI8apw/wzKfb3I7InE1cFMztx0kNpFfcBAZ1asXIy5u6HVW+iWrZOEfq8Xg0IiLC7TCMKRzxB+G1MHaVv5Sro8fy3ugOtCtBp0DKjMQ/4e2upPwZxQ0nH+ai1u2ZMiAEn2L6LAeAiKxXVU/W9faEuTGlQZV6cPXDNIv7gZur/MrDSzaRnGqFE4uVlNOw4CbSj+1kZOK9BJ0fyvP9WhfrxHE2ljyMKS3CboM6rXjYdw6HjhzhzW93uR2ROSM9DT66FfZ9z/2pd/BncAemDW1LgF/J/RFcciM3xvyVUzgxMDGG/wav4JW1O6xwYnGgCisnwdaPeUmG81PFf/DOyDAqBbo3haw3WPIwpjRp0BY8o7g6fgktZQ//XraZsnJds9j69mX4eQbv+d3IXOnBnFHh1K5czu2oCsyShzGlzTWPIBVqMa3aXL7cdphVm6Pdjqjs2rAA1j7GV4FX8mjSIN4aEUazoEpuR+UVljyMKW3KV4Pr/lc48bGPrXCiKyLXosvGsjkwhNvjR/HqYA+hjbItPF4iWfIwpjS6tD807cSdae+SEhfNfz77w+2IypaDG9BFwzjo34SBcWN5tHconVsEux2VV1nyMKY0OlM4MS2JmXUW8/Z3e9hy0AonFonju2FeP+KlMr3jxnNr59YMCm/kdlReZ8nDmNKqVnO4fBxtYlfTpdxWHlqy0QonFraTR2FuH5KSk+kTP4Frwltx7zXN3Y6qUFjyMKY06zgRqjflhfJz2LwvhoXrrHBioUk+CfMHkBZ7gJtOjqfpRW14ouclxXIKWW+w5GFMaeZfDrq/SOWTe3ii1mc8t3IbR61wovelpcL7I9GDv3J3yt2kNwjn1cGh+PmW3h+xpXdkxpgM53eGln3on7iIWslRPL1iq9sRlS6qsHwc7FjFk3or26p25K3hYZQPKL5TyHqDJQ9jyoLrnsbHL5A3ay3ko1+i+GHnMbcjKj2+fBZ+fZd3fPuzzP86Zo8Kp0bFALejKnSWPIwpC6rUhasn0zTuJ4ZX+YWHl2y0woneEPEOfPUsqwI682JKX94ZEUbDGhXcjqpIWPIwpqwIuwXqhvAvn3c5EhPDzG+scGKBbFuBfjKBXwPDuPfEcKbf7OGS+lXdjqrIFCh5iEgNEVkjIjuc92wfnxSRriKyXUQiRWRSbv1FpImIJGaaCGp6pj5tRWSjs69XpLTeymCMt/n4Qo+XCTwdwyvBn/DK2h3sO2aFE8/J/p/RD0axN+ACboq7g2f6h9KxeZDbURWpgn7zmASsVdXmwFpn+S9ExBd4HegGtAAGi0iLPPTfqaohzmtMpvXTgNFAc+fVtYBjMKbsqB8KYbdyVfxSWvns5pFlm6xwYn4d3QHzB/Knb036xt3Lvd1C6N2mgdtRFbmCJo+ewGzn82ygVzZtwoFIVd2lqsnAQqdfXvv/PxGpC1RR1R8041/8nNz6GGOyuGYyUjGIqVXf5evt0azcdNjtiEqOhMPwbh8SU6FX/ERuuKw1t3dq5nZUriho8ghW1UMAznvtbNrUBzI/mRTlrMutf1MR+VVEvhKRjpn2FZXDvv5GREaLSISIRMTExORnXMaUXuWqwnVPE5Swhfuqf8NjH2/hhBVOzN3peJjbj9QTRxlwYgKXXhLCIz1alNqHAHOTa/IQkc9EZFM2r5659T2zi2zW5fY9+RDQSFXbABOA+SJSJb/7UtUZqupRVU9QUNk6H2nMWV3SF5pdxe2p80lPOMTLa6xw4lmlJsN7Q0k/spXRSXdToYmHlwaU3ClkvSHXqaxUtXNO20QkWkTqquoh55TSkWyaRQENMy03AA46n7Ptr6pJQJLzeb2I7AQucPbVIId9GWPyyimc6Du1AzODF9P7u+r0blO/TN0tlGfp6bD0Ttj9FQ/rnUTVvJz3h3ko51+6HwLMTUFPWy0DhjufhwNLs2mzDmguIk1FJAAY5PTLsb+IBDkX2hGRZmRcGN/lnNpKEJH2zl1Ww3I4pjEmNzXPg44TaB37Gd3Kb+XhJZuscGJ2PnsENr7PVJ+b+KJcZ2aPCqdqeX+3o3JdQZPHs0AXEdkBdHGWEZF6IrICQFVTgbHAKmArsEhVN5+tP9AJ+F1EfgM+AMao6nFn2x3Am0AksBP4tIBjMKbsunwc1GjGc+XnsHX/ERas2+d2RMXLD1Ph+1dZ4t+N6ek9mTUynLpVy7sdVbEgZeU2PY/HoxEREW6HYUzxs/NzeLc3H1QayuMnbmTtxKsIqhzodlTu2/QR+sEofgrswIiTY5k9qj3tmtV0O6oiJyLrVdWTdb09YW5MWXfe1XBJP/qeWkTtlCiescKJsPsbdPHtRAa2ZET8aF4eGFomE8fZWPIwxsB1TyP+5Xiz1nt89GsU3+886nZE7onejC68iRi/evSLu4dJPULodmldt6Mqdix5GGOgcjBc8whN4n5iRJVfeHjJJpJS09yOqujFRcHcfpzUcvSKm8DgK1sz4vKmbkdVLFnyMMZk8IyCem140GcOR2OOMPPrMlY48dRxmNuX5MR4+iZMoH2b1jzQ9UK3oyq2LHkYYzKcKZyYdIxXgj/h1c8jy07hxJREWHgT6cd2MiJxPLXPD+W5fq3K7NPjeWHJwxjzP/XaQNhtXBm3lNY+u5i8tAwUTkxPg49uQ/f9yMTUu4iv045pQ9viX4qnkPUG+9MxxvzV1Q8hlWoztepcvvkjmk9Lc+FEVfj0Adj6MS/JcCIqXcnbI8KoFJhr8Y0yz5KHMeavylWFrs9QK2EL99f4hsc+3lx6Cyd++zKsm8kC/17M9+nBnFHtqF25nNtRlQiWPIwxf9eyDzT7B6NT5kHCYaasLoWFEzcsgLWP8WXgVTxxeiBvDffQtFZFt6MqMSx5GGP+TgS6v4Rvegozgz9i1ve72XQgzu2ovCfyM3TZWDYHtmFMwi28NqQtbRplOxGqyYElD2NM9mqeBx0n0ip2Ld3Kb+GhJZtIKw2FEw/+ir43jIP+TRgYdxeP9Q7h6ouC3Y6qxLHkYYzJ2RXjoOb5PFd+Ntv2H2H+zyW8cOLxXTCvP/E+VegVN57bOocwMKyR21GVSJY8jDE58wuE7i9R6eR+nglaw/MrtxGTkOR2VOfm5FGY25ek5GR6x0+kc3gr7rnmfLejKrEseRhjzq7ZVXBpf3qffJ+6KVE89ckWtyPKv+STMK8/abEHGHxyPM0uCuGJnpfYQ4AFYMnDGJO7a59C/Mszs9YClmw4wHeRJahwYloKvD8CPbSBu1LuhgbhvDo4FD97CLBA7E/PGJO7ysHQ+REax63jlqoRTC4phRNVYfk42LGaJ/QW/qjWkbeGh1E+oGxPIesNBUoeIlJDRNaIyA7nPdt73USkq4hsF5FIEZmUW38RGSIiGzK90kUkxNn2pbOvM9tqF2QMxpg8ajsS6oXygMzh6NEjTP+yBBRO/OJp+HUub/v252P/rsweGU71igFuR1UqFPSbxyRgrao2B9Y6y3/hzEX+OtANaAEMFpEWZ+uvqvNUNURVQ4CbgT2quiHTboec2a6qRwo4BmNMXjiFEwOS/uS14OW8/mUke46edDuqnEW8DV8/z0r/LkxJ7ceskWE0rFHB7ahKjYImj57AbOfzbKBXNm3CgUhV3aWqycBCp19e+w8GFhQwTmOMN9QLgfDb6Ri3jFDfYlw4cety9JOJ/BIYxvhTw3jjZg8t61V1O6pSpaDJI1hVDwE479mdQqoP7M+0HOWsy2v/gfw9ebzjnLKaLGe5XUJERotIhIhExMTE5G1Expiz+8e/kMp1eL3KHL7fEc0nGw+5HdFf7fsJ/fAW9gZewJC4O3i2f1suP7+W21GVOrkmDxH5TEQ2ZfPqmVvfM7vIZl2eflURkXbAKVXdlGn1EFW9FOjovG7Oqb+qzlBVj6p6goKC8hiuMeasylWBrs9QM2EbD9T4hsc/3kL86RS3o8oQ8we6YCDHfWvRJ3Yc468PoWdI/dz7mXzLNXmoamdVvSSb11IgWkTqAjjv2V1/iAIaZlpuABx0PufWfxBZvnWo6gHnPQGYT8ZpMWNMUWrRC87vzC0p8/E5cah4FE6MPwRz+5CYKvSKn0jPy1txW8dmbkdVahX0tNUyYLjzeTiwNJs264DmItJURALISAjLcusvIj5AfzKukZxZ5ycitZzP/kAPIPO3EmNMURCB61/AV1N5M/gj5vywh41RLhZOPB0P8/qTevIYA05MoNUlIUzu3sIeAixEBU0ezwJdRGQH0MVZRkTqicgKAFVNBcYCq4CtwCJV3Xy2/o5OQJSqZr4fMBBYJSK/AxuAA8DMAo7BGHMuajSDjvdxSezn3FBhMw8t2ehO4cTUZHhvCOlHtnLb6Xup2MTDSwNa4+NjiaMwSbG8U6IQeDwejYiIcDsMY0qX1CSYdjknE08TevwJHu7Zhps7NCm646enw0e3waYPeFDH8mv163jv9g5ULe9fdDGUciKyXlU9WdfbE+bGmHPnFE6seGo/zwWt5vmV2zmScLrojr9mMmz6gNd8hvBVuauZNTLcEkcRseRhjCmYZldCq4H0PPk+9dP28+TyrUVz3B+mwg+vsdi/OzPSb2T2qHDqVLUpZIuKJQ9jTMFd+yQSUIEZNRew7LcDfLujkAsnbvoQVj3ID4GX82DiEN4aEU7z4MqFe0zzF5Y8jDEFV6k2XPNvGsVFcFvVdUxeuonTKYVUOHH31+jiMfxR7lJGxN/GfwaFEtakRuEcy+TIkocxxjvajoT6Hu6Xdzl2NJrpX+30/jEOb0IXDuGIXz36xd7NQze2oesldb1/HJMrSx7GGO/w8YEeUwhI+pPXg5cz9cud7PZm4cTY/TCvHye1HL3iJjL0qtYMK8o7u8xfWPIwxnhP3dbQbgxXxH2Mx3cnk5d4qXDiqeMwty/JiSfomzCBDqGtuP+6Cwu+X3POLHkYY7zrH/9CKtfl1Spz+CEymo9/L2DhxJREWDCI9OO7GJE4juDmbXmubyt7etxlljyMMd4VWNkpnLidSTW+5onlBSicmJ4GH96K7v+Z8Sl3kVCnPdOGhOJvU8i6zv4GjDHe16InnN+FUSnz8T1xiJdWbc//PlRhxf2wbTkvygh+rXwVb48Io2Kgn/fjNflmycMY433/XzgxjbeCP2TOj3v5PSo2f/v45iWIeIv5fr1Y4NOd2aPCCaocWCjhmvyz5GGMKRw1mkKn+2gZ+wU3VtjMQ4s35b1w4q/z4PMn+CLgKp5MGsjbI8JoWqti4cZr8sWShzGm8Fx2D9S6gKcDZ/PHgRjm/rg39z471qDL7mZTYBvuOHELrw1pS0jDaoUeqskfSx7GmMLjFwjdp1DxVBTP117Fi6u2cyT+LIUTD6xHFw3nQEBTBsXdxeO923D1RcFFF6/JM0sexpjC1bQjtBrEjSc+oH7afp74JIfCicd2wrwBxPlUoXfceEZ3CWFAWMPs2xrXWfIwxhQ+p3DizJrz+fi3A3z9R8xft5+Igbl9OZ2SSp/4iXRp15q7rz7fnVhNnhQoeYhIDRFZIyI7nPfqObTrKiLbRSRSRCZlWt9fRDaLSLqIeLL0edBpv11Ersu0vq2IbHS2vSL2pJAxxV+lIOj8GA3j1nN7tXXcNe8Xvth2JGNb0gl0fn9S4w5x08lxNLuoDY/f2NIeAizmCvrNYxKwVlWbA2ud5b8QEV/gdaAb0AIYLCItnM2bgD7A11n6tCBjrvOWQFdgqrMfgGnAaKC58+pawDEYY4pC6HBoEMY/5V1aVE9n1Ox1TFu7lbg5Q0g/+Bu3n74Lv8bteXVwG/zsIcBir6BP2/QErnI+zwa+BB7I0iYciDwzF7mILHT6bVHVrc667Pa7UFWTgN0iEgmEi8geoIqq/uD0mwP0Aj4t4DiMMYXNxwd6vIzvG1ey0G88hysHkvb1CarKUZ7yvYPreo+kb9sG+Nrc4yVCQZNHsKoeAlDVQyJSO5s29YH9mZajgHa57Lc+8GOWPvWBFOdz1vXZEpHRZHxLoVGjRrkc0hhT6OpcCr2mIds/oQ6w99gpPq/SjnH9xtuT4yVMrn9bIvIZUCebTQ/l8RjZ/RqR25NCOfXJ175UdQYwA8Dj8XihtKcxpsBaD4TWAxGgifMyJU+uyUNVO+e0TUSiRaSu862jLnAkm2ZRQOb77RoAB3M5bE59opzP+dmXMcYYLyvoVallwHDn83BgaTZt1gHNRaSpiASQcSF8WR72O0hEAkWkKRkXxn92TpEliEh75y6rYTkc0xhjTCEqaPJ4FugiIjuALs4yIlJPRFYAqGoqMBZYBWwFFqnqZqddbxGJAjoAn4jIKqfPZmARsAVYCdylqmcmRL4DeBOIBHZiF8uNMabIiVdm+SoBPB6PRkREuB2GMcaUKCKyXlU9WdfbzdTGGGPyzZKHMcaYfLPkYYwxJt8seRhjjMm3MnPBXERigDzMRJOtWsBRL4ZTEtiYy4ayNuayNl4o+Jgbq2pQ1pVlJnkUhIhEZHe3QWlmYy4bytqYy9p4ofDGbKetjDHG5JslD2OMMflmySNvZrgdgAtszGVDWRtzWRsvFNKY7ZqHMcaYfLNvHsYYY/LNkocxxph8s+RxFiLSVUS2i0ikiPxtfvbSQEQaisgXIrJVRDaLyL3O+hoiskZEdjjv1d2O1dtExFdEfhWR5c5yqR6ziFQTkQ9EZJvz992hDIx5vPPvepOILBCRcqVtzCLytogcEZFNmdblOEYRedD5mbZdRK471+Na8siBiPgCrwPdgBbAYBFp4W5UhSIVmKiqFwPtgbuccU4C1qpqc2Cts1za3EvGNAFnlPYx/xdYqaoXAa3JGHupHbOI1AfuATyqegngS8Z8QqVtzLOArlnWZTtG5//2IKCl02eq87Mu3yx55CwciFTVXaqaDCwEerock9ep6iFV/cX5nEDGD5T6ZIx1ttNsNtDLlQALiYg0ALqTMTfMGaV2zCJSBegEvAWgqsmqGkspHrPDDygvIn5ABTJmHi1VY1bVr4HjWVbnNMaewEJVTVLV3WTMixR+Lse15JGz+sD+TMtRzrpSS0SaAG2An4BgZ+ZGnPfaLoZWGP4D/BNIz7SuNI+5GRADvOOcqntTRCpSisesqgeAF4F9wCEgTlVXU4rHnElOY/TazzVLHjmTbNaV2vuaRaQS8CEwTlXj3Y6nMIlID+CIqq53O5Yi5AeEAtNUtQ1wkpJ/uuasnPP8PYGmQD2googMdTcq13nt55olj5xFAQ0zLTcg4ytvqSMi/mQkjnmq+pGzOlpE6jrb6wJH3IqvEFwO3Cgie8g4HXm1iMyldI85CohS1Z+c5Q/ISCalecydgd2qGqOqKcBHwGWU7jGfkdMYvfZzzZJHztYBzUWkqYgEkHGRaZnLMXmdiAgZ58G3quqUTJuWAcOdz8OBpUUdW2FR1QdVtYGqNiHj7/VzVR1K6R7zYWC/iFzorLoG2EIpHjMZp6vai0gF59/5NWRc0yvNYz4jpzEuAwaJSKCINAWaAz+fywHsCfOzEJHryTg37gu8rapPuRuR94nIFcA3wEb+d/7/X2Rc91gENCLjP2F/Vc16Ua7EE5GrgPtUtYeI1KQUj1lEQsi4QSAA2AWMJOMXyNI85seAgWTcVfgrcCtQiVI0ZhFZAFxFRun1aODfwBJyGKOIPASMIuPPZJyqfnpOx7XkYYwxJr/stJUxxph8s+RhjDEm3yx5GGOMyTdLHsYYY/LNkocxxph8s+RhjDEm3yx5GGOMybf/A/G9gQTd4LvUAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.plot(np.fft.ifft(1j * y0_h.imag))\n", "plt.plot(y0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our equation becomes\n", "\n", "$$\\frac{\\partial \\hat{y}}{\\partial t} = \\hat{z} $$\n", "\n", "$$\\frac{\\partial \\hat{z}}{\\partial t} = -c^2 \\gamma \\hat{z} - c^2 k ^2 (1+k^2 l^2)\\hat{y} $$\n", "\n", "Let $S=(\\hat{y}, \\hat{z})$" ] }, { "cell_type": "code", "execution_count": 144, "metadata": {}, "outputs": [], "source": [ "def dSdt(S, t, k):\n", " y_h_ri = S[:2*N]; z_h_ri = S[2*N:]\n", " y_h = 0*y_h_ri[:N] + (1j) * y_h_ri[N:]\n", " z_h = 0*z_h_ri[:N] + (1j) * z_h_ri[N:]\n", " dy_h = z_h\n", " dz_h = -c**2 * gamma * z_h - c**2 * k**2 * (1+(k**2) * (l**2)) * y_h\n", " dy_hri = np.concatenate([dy_h.real, dy_h.imag])\n", " dz_hri = np.concatenate([dz_h.real, dz_h.imag])\n", " return np.concatenate([dy_hri, dz_hri])" ] }, { "cell_type": "code", "execution_count": 145, "metadata": {}, "outputs": [], "source": [ "dt = 5e-6\n", "t = np.arange(0, 0.025, dt)\n", "sol = odeint(dSdt, y0=S0, t=t, args=(kappa,))" ] }, { "cell_type": "code", "execution_count": 146, "metadata": {}, "outputs": [], "source": [ "y_h = sol.T[:N] + (1j) * sol.T[N:2*N]\n", "y = np.zeros_like(y_h)\n", "for k in range(len(t)):\n", " y[:,k] = np.fft.ifft(y_h[:,k]).real\n", "y = y.real" ] }, { "cell_type": "code", "execution_count": 147, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[-6.50521303e-19, -2.16840434e-18, -4.33680869e-19, ...,\n", " -8.67361738e-19, -4.33680869e-19, -1.08420217e-18],\n", " [ 3.44827586e-04, 3.50354508e-04, 3.66091098e-04, ...,\n", " 5.17285057e-04, 5.17288693e-04, 5.17270589e-04],\n", " [ 6.89655172e-04, 6.88798949e-04, 6.86874136e-04, ...,\n", " 8.61976570e-04, 8.61969122e-04, 8.62007466e-04],\n", " ...,\n", " [-6.89655172e-04, -6.89927843e-04, -6.90511478e-04, ...,\n", " -8.62221702e-04, -8.62233294e-04, -8.62169773e-04],\n", " [-3.44827586e-04, -3.43971363e-04, -3.42046549e-04, ...,\n", " -5.17148984e-04, -5.17141536e-04, -5.17179880e-04],\n", " [ 2.87151175e-18, -5.52692199e-06, -2.12635117e-05, ...,\n", " -1.72457471e-04, -1.72461107e-04, -1.72443003e-04]])" ] }, "execution_count": 147, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y" ] }, { "cell_type": "code", "execution_count": 151, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAD4CAYAAAAZ1BptAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAvoklEQVR4nO3deXxV1b3//9fnnJN5JIQkhDCFQSaZEpNQJ6hapdpibZVBAVEIXvXett/e69DpDm3vt7e13tbaKgFkUCFQq1e+yq21aKxVEkgCMo8hCUMYM0ASQqb1+yMHfzFNQpIz7JOTz/Px4JGcc/ba+7MY8mbvtfdaYoxBKaWUao/N6gKUUkr5Lg0JpZRSHdKQUEop1SENCaWUUh3SkFBKKdUhh9UFuFNsbKwZNmxYj9vX1NQQFhbmvoJ8XF/rL2if+wrtc/cUFBScN8YMaO8zvwqJYcOGkZ+f3+P2OTk5TJ8+3X0F+bi+1l/QPvcV2ufuEZGSjj7Ty01KKaU6pCGhlFKqQxoSSimlOqQhoZRSqkMaEkoppTqkIaGUUqpDGhJKKaU6pCHhJR8ePMu7u8qoa2iyuhSllOoyv3qYzledu3SFpa8WUN/YTGSwg69NSuTRm4aTPCDc6tKUUqpTeibhBWs+LaahqZnn7p/El8fE8cfCE3zzpU85fOaS1aUppVSnNCQ8rOZKI2u3FnPnuAS+lZLEr+dM4b3v3EKA3caDK/IovVBrdYlKKdUhDQkPy95+nIt1jSy9Nfnz94b2D+O1xenUNzUzb0Uup6vqLKxQKaU6piHhQQ1Nzaz8uIi04TFMGdLvC5+Njo9g7SNpVNY2MCdrKycrL1tUpVJKdUxDwoPe2XWKU1V1PNbqLKK1iUnRrHkkjQs19Tzw8laKz9d4uUKllOqchoQHrfzbMUbHhzN9dFyH26QM7cf6JRnU1jdy/7KtHNLBbKWUD9GQ8JDz1VfYc/Ii901NwmaTTredMCiKDUunATAnK5cDpy96o0SllLomDQkP2X6sHID04TFd2n50fAQbl04jwC7MW57H/jINCqWU9TQkPCTvWDkhAXYmDIrqcpvhsWFkZ04j0G5j3vJc9p3SoFBKWUtDwkPyjpWTMrQfAfbu/Ra3BEUGQQ4781bksudklYcqVEqpa9OQ8ICq2gYOnL5IWhcvNbU1LDaMDUszCAt0MG95LjuPV7q3QKWU6iINCQ/ILynHGHocEtDywN2GpRlEhwby0Io88ovL3VihUkp1jYaEB2w7Vk6g3cbkwdEu7SepXygbl04jLiKIBa9sI6/ognsKVEqpLtKQ8IC8Y+VMGhxFcIDd5X0lRAWTvTSDgVHBPLxqO1uPalAopbxHQ8LNaq40svtklUuXmtqKiwgmO3MaSf1CWLR6G58eOe+2fSulVGc0JNyssLSCpmZD+vD+bt3vgIgg1mdmMCQmlEWrt/PXQ+fcun+llGqPhoSbbTtWjt0mTB3a79obd1NseBDrl2SQPCCcxWvy+eDAGbcfQymlWtOQcLO8Y+VMSIwkPMgzi/71Dw9i/ZJ0rkuIYOmrBfxpz2mPHEcppUBDwq2uNDax83glNwxz33hEe6JDA3ltcToTBkXxxLpCNu8u8+jxlFJ9l4aEG+05eZH6xmZSh7n/UlNbUSEBvPpoOlMGR/OP63fwzq5THj+mUqrv0ZBwo8KSCgCPjEe0JzzIwepH0pg6JJpvZ+9k02caFEop99KQcKOCkgqGxIQSFxHstWOGBzlYvSiNlCH9+E72Dt4sPOG1Yyul/J+GhJsYY8gvqSDFS2cRrYUFOVj9yA1kJPfne3/4jOxtpV6vQSnlnzQk3OR4+WXOV1+xJCQAQgMdvPLwDdwyagDPvLmbtVuLLalDKeVfNCTcpKC0ZQI+q0ICIDjATtaCFG4fG8+P397Lqk+OWVaLUso/uCUkROQuETkoIkdE5Jl2PhcRecH5+S4RmXqttiLySxE54Nz+LRGJdketnpJfXEFEkIPR8RGW1hHksPP7B6dy5/h4/v3/7eOVv2lQKKV6zuWQEBE78DtgJjAOmCsi49psNhMY5fyVCbzUhbbvAxOMMROBQ8CzrtbqSQUlFUweEo39GutZe0Ogw8aL81qC4j/e2cdKDQqlVA+540wiDThijCkyxtQD2cCsNtvMAtaaFrlAtIgM7KytMebPxphGZ/tcIMkNtXrEpboGDp65ZOmlprYC7C1BMXNCAj95Zx8vf3TU6pKUUr2QO+aOGAQcb/X6BJDehW0GdbEtwCPAhvYOLiKZtJydEB8fT05OTjdK/6Lq6uoetd9zvhFjIKDyODk5vvWswjcTDeXn7fz8fw+w/9BR7h0ZgEjL2U5P+9ubaZ/7Bu2z+7gjJNq7vmK6uM0124rID4BG4PX2Dm6MyQKyAFJTU8306dOvUW7HcnJy6En7He8fwiaHWXDPLUQEB/T4+J4yY7rh2Td3sTH/BPGDBvPszDGISI/725tpn/sG7bP7uCMkTgCDW71OAtr+d7qjbQI7aysiC4F7gNuMMW2Dx2cUllZwXUKkTwYEgN0m/Py+iYQE2Mn6axHGGL7/1bFWl6WU6gXcERLbgVEiMhw4CcwB5rXZZhPwpIhk03I5qcoYUyYi5zpqKyJ3AU8Dtxpjat1Qp0c0NRt2lFZy75REq0vplM0m/NvXxwOw/ONjOOw20oJ8NneVUj7C5ZAwxjSKyJPAe4AdeMUYs1dEHnN+/jKwGfgqcASoBRZ11ta56xeBIOB95zX0XGPMY67W624HTl+k+kojqUM9O/OrO4i0BEVDs+GlnKOcGhHA9Onm8zEKpZRqyy2LHhhjNtMSBK3fe7nV9wZ4oqttne+PdEdtnnZ1Uj9furOpMyLCT2dNoLGpmY35J4h5Zx8/vHucT9y6q5TyPfrEtYvySyqIjwwiqV+I1aV0mc05RnHHUAerPinmyXWF1DU0WV2WUsoHaUi4KL+4gtShMb3uko3NJjw4Nogf3j2W/91zmodW5FFzpfHaDZVSfYqGhAtOV9VxsvKy19aP8ITFNyfz27lTKCit4JfvHbS6HKWUj9GQcEF+Scukfqm9OCQAvjYpkYXThrH602K2HSu3uhyllA/RkHBBQUkFIQF2xiVGWl2Ky5666zoGx4Tw1BufcblexyeUUi00JFxQUFLBpMFRBNh7/29jaKCD/7pvIsUXann+fb3spJRq0ft/ulmktr6Rvacu9ppbX7viSyNjmZc+hJV/O8a+UxetLkcp5QM0JHpo5/FKmppNr3iIrjuevmsMQQ47qz/V6cWVUhoSPXb1IbqpQ/znTAIgKiSAb0wdxNs7T1FRU291OUopi2lI9FB+SQWj48OJCvXNSf1csWDaUK40NrMh//i1N1ZK+TUNiR5objYUllT41XhEa2MSIslIjuHVrSU0NeskgEr1ZRoSPXD4bDUX6xpJ8bPxiNYWThvGycrLfHDgrNWlKKUspCHRAwXO8Yje/hBdZ+4YF8/AqGDWbi22uhSllIU0JHogv6Sc/mGBDO0fanUpHuOw23gwfQgfHz7P/jK9HVapvkpDogcKnOMRvW1Sv+6alz6UmLBAvrthp84Sq1QfpSHRTecuXaHkQi2pw/z3UtNVMWGB/Or+SRw4fYn/u3m/1eUopSygIdFNBZ8vMuS/g9atzRgTx6M3DWfN1hLe33fG6nKUUl6mIdFNBSXlBDpsTBjU+yf166qn7rqO8YmR/Msbn1F6wWeXG1dKeYCGRDfll1QwcVAUQQ671aV4TZDDzm/nTsEY+ObLn7L3VJXVJSmlvERDohvqGprYc7KKlD4wHtFW8oBw3nhsGgE2YfayXD49ct7qkpRSXqAh0Q27T1bR0GRI8bP5mrpqVHwEf3z8SwyKDmHhqm38ac9pq0tSSnmYhkQ35BdfHbTumyEBMDAqhI2PTeP6QVE8ua5QB7OV8nMaEt1QUFJOcmwY/cODrC7FUlEhAax+JI3xg6J4/PUCPjigQaGUv9KQ6CJjDAUlFUztw2cRrUUGB7D2kTSuS4jgsVcL+fCgzvGklD/SkOiiovM1VNQ2+PV8Td0VFRLAa4+mMyo+nMy1+TpGoZQf0pDoogLneERfeNK6O6JDA1m3JIMJg6J4Yl0hb+88aXVJSik30pDoovyScqJCAkiODbe6FJ8TFRLAq4+mkzq0H9/ZsJM/6GJFSvkNDYkuyndO6mez+fekfj0VHuRg9aI0bhoZy1N/3MXG7RoUSvkDDYkuKK+pp+hcTZ++9bUrQgLtLF+Q+nlQbNheanVJSikXuSUkROQuETkoIkdE5Jl2PhcRecH5+S4RmXqttiJyv4jsFZFmEUl1R509VdgHFhlyl+CAlqC4dfQAnv7jbtZv06BQqjdzOSRExA78DpgJjAPmisi4NpvNBEY5f2UCL3Wh7R7gPuCvrtboqvySChw2YWJStNWl9ArBAXaWzU9hxnUDePbN3bq6nVK9mDvOJNKAI8aYImNMPZANzGqzzSxgrWmRC0SLyMDO2hpj9htjDrqhPpcVllQwflAUIYF9Z1I/VwUH2Hl5fgp3jIvnx2/vZcXHRVaXpJTqAYcb9jEIaD1KeQJI78I2g7rYtlMikknL2Qnx8fHk5OR0p/kXVFdX/137xmbDjtJavjzY4dK+fVF7/XW32UmGynI7P313P2UlR7k5KcCjx7sWb/TZ12if+wZP9dkdIdHe7T6mi9t0pW2njDFZQBZAamqqmT59eneaf0FOTg5t2xeWVtDw50+ZddNEpl8/sMf79kXt9dcTpt/azNde/IT8Khs/euhGjx+vM97qsy/RPvcNnuqzOy43nQAGt3qdBJzq4jZdaWspHbR2ncNu456JA/nseCWnq+qsLkcp1Q3uCIntwCgRGS4igcAcYFObbTYBC5x3OWUAVcaYsi62tVR+cQWDY0KIiwy2upRe7c7x8QC8v0+n7lCqN3E5JIwxjcCTwHvAfmCjMWaviDwmIo85N9sMFAFHgOXA4521BRCRb4jICWAa8K6IvOdqrd1ljCG/pILUPrKetSeNGBBOcmwY7+3VGWOV6k3cMSaBMWYzLUHQ+r2XW31vgCe62tb5/lvAW+6or6eOl1/mfPUVnfnVDUSEr4xPYMXHRVTVNhAVau0AtlKqa/SJ607kl5QDOh7hLneOj6ex2bBF159QqtfQkOhEfkkFEUEORsdHWF2KX5iUFE18ZBB/1ktOSvUaGhKdKCypYPKQaOw6qZ9b2GzCV8Yl8NGhc9Q1NFldjlKqCzQkOlB1uYGDZy7poLWb3Tk+gcsNTfz10DmrS1FKdYGGRAd2lFZgjC4y5G7pyTFEBjvYsl+XO1WqN9CQ6EBhSQU2gcmDo60uxa8E2G2kDe/PtuJyq0tRSnWBhkQH8ksqGDswkrAgt9wlrFpJHx7DsfM1nL2oT18r5es0JNrR2NTMzuOVeuurh6QNbxnn0bMJpXyfhkQ7Dpy+RG19EynDdNDaE8YnRhIaaGfbMQ0JpXydhkQ78ov1ITpPcthtpAztpyGhVC+gIdGO/JIKBkYFkxgdYnUpfit9eAwHTl+isrbe6lKUUp3QkGhHYUkFKXoW4VFpw/sDsL24wuJKlFKd0ZBoo6zqMqeq6jQkPGxiUhSBDhvbjl2wuhSlVCc0JNo4ePoSAOMToyyuxL8FB9iZnBTNNj2TUMqnaUi0UXSuBoDkAWEWV+L/0obHsOdkFTVXGq0uRSnVAQ2JNorOVxMZ7KB/WKDVpfi9tOExNDUbCkv1bEIpX6Uh0UbRuRqSB4QjojO/etrUof2w24S/HT5vdSlKqQ5oSLTREhJ6qckbwoMcfHlMHMs/LuLtnSetLsdSxhi+k72DZ9/cpdOVKJ+iIdHK5UbD6Yt1jBgQbnUpfcZv5kwmbXgM392wkzcLT1hdjmXe3V3G/+w8Rfb240x/LocXPzisa24on6Ah0cqZmmYAkmP1TMJbQgMdrHo4jWkj+vO9P3zGxvzjVpfkdfWNzfzyvYOMSYhgy/+5lZtHxfLcnw8xJyuXah3UVxbTkGilrMYAMCJOzyS8KSTQzsqFN3DTyFieemMX6/JKrS7Jq9bllVByoZanZ44heUA4y+an8tKDU9l9sopHVm/ncr2eUSjraEi0crqmGZvA0P6hVpfS5wQH2Fm+IJUZ1w3g+2/tZu3WYqtL8opLdQ288MERvjSiP9NHD/j8/ZnXD+S/Z09me3E5ma/mc6VRg0JZQ0OilbKaZpL6hRLksFtdSp8UHGDn5fkp3D42nh+/vZeVfztmdUket+yjIspr6nl25ti/u6Pu65MS+a9vTuTjw+d5/LVCHaNQltCQaOV0jdE7mywW5LDz+wenMnNCAj95Zx9Zfz1qdUke09DUzOpPi7l74kCuT2r/Cf8HUgfzs29MYMuBsyxek09tvY5RKO/SkHBqbjacrm0mOVbHI6wW6LDxwtwp3D1xIP+5+QC/zzlidUkesetEFdVXGrnn+oGdbvdg+lCeu38Snx49z8JXtnGprsFLFSqlIfG5sot11DfpdBy+IsBu4zezJzNrciK/+NNBfrvlsNUluV1uUcvkhldX6uvMt1KSeGHuFHaUVvLQym1UXdagUN6hIeFUdK4a0JDwJQ67jecfmMx9Uwfxq/cP8fz7hzDGWF2W2+QWXeC6+Aj6hwd1aft7Jiby+wensu9UFQtW5mlQKK/QkHC6OrHfSH2QzqfYbcIvvzWJB1KTeGHLYX753kG/CIqGpmbyiyvISO7eErlfGZ/ASw+msK/sIvNX5lFVq0GhPEtDwqnoXDXBdhgQ0bX/1SnvsduEn983kXnpQ/h9zlH+c/P+Xh8Uu05UcrmhiYzk/t1ue/u4eF5+KIUDZZd4cGUuFTW6up/yHA0Jp6LzNQwMs+nEfj7KZhN+du8EHv7SMJZ/fIz/eGdfrw6K3KKW9b3TexASALeNjWfZghQOnalm7vJcLlRfcWd5Sn3OLSEhIneJyEEROSIiz7TzuYjIC87Pd4nI1Gu1FZEYEXlfRA47v3p0qbijZ6tJCNOA8GUiwr9+bRyP3jScVZ8U86+b9tLc3DuDIrfoAmMSIohxYUr6GdfFsXJhKsUXapiTlcvZSzoxoHI/l0NCROzA74CZwDhgroiMa7PZTGCU81cm8FIX2j4DbDHGjAK2OF97RG19I6eq6kgI0xMrXyci/PDusSy9JZm1W0v4wf/s6XVBUd94dTyiZ2cRrd08agCrHk7jZOVl5mTlckZnkFVu5o6fimnAEWNMkTGmHsgGZrXZZhaw1rTIBaJFZOA12s4C1ji/XwPc64Za23XsfMug9UANiV5BRHhm5hgenz6C9dtKeebNXTT1oqDYffLqeET3Bq07Mm1Ef9Y+ksaZqjrmZOVSVnXZLftVCsDhhn0MAlpP3XkCSO/CNoOu0TbeGFMGYIwpE5G49g4uIpm0nJ0QHx9PTk5OtztQXtfMt0YHkBBQ16P2vVV1dXWv7u8NQYZZIwLYmH+CE6dOs/j6QGzXGFPyhT5vOtoy0Nxw6gA55w+6bb/fnRLArwpq+PqvP+SZtGD6h7T8p8cX+uxt2mf3cUdItPevsu1/6zrapittO2WMyQKyAFJTU8306dO70/xz9wE5OTn0tH1v5A/9nTEDRmw5zPPvHyJ2QBzPPzAJh73jM0Jf6POKI3mMSbjC175yi1v3Ox1ITa1k/so8/nsXrF+SxuCYUJ/os7dpn93HHddXTgCDW71OAk51cZvO2p5xXpLC+fWsG2pVfuifbhvF03eNYdNnp/h29k4ampqtLqlDDU3N5JeUu2U8oj2TB0ezbnEGl+oamb1sKyUXajxyHNV3uCMktgOjRGS4iAQCc4BNbbbZBCxw3uWUAVQ5LyV11nYTsND5/ULgbTfUqvzUP0wfwQ/vHsu7u8t4cl0h9Y2+GRT7yy5S19BM6jDP3ax3fVIU65akc7mhidnLcjld45u/F6p3cDkkjDGNwJPAe8B+YKMxZq+IPCYijzk32wwUAUeA5cDjnbV1tvk5cIeIHAbucL5WqkOLb07mX782jvf2nuHx1wt8cg2GgpIKAFKGevSObsYnRrE+M4OGpmb+77Y6jpy95NHjKf/lltt5jDGbjTGjjTEjjDE/c773sjHmZef3xhjzhPPz640x+Z21db5/wRhzmzFmlPNruTtqVf5t0Y3D+cm9E/jL/rM89mqBz63BUFBSQWJUMAOjQjx+rDEJkWRnZmAMzMnK5dAZDQp/9fbOkxws98zfdb3nU/md+RlD+c9vXM+HB8+xZG2+TwVFYUkFUz18FtHaqPgInkkLxibCnKxc9pdd9NqxlXcYY/jpu/vJOe6Zebw0JJRfmpc+hF98ayJ/O3KeR9f4xjrRpyovc6qqzuOXmtpKDLexYek0ghw25i7PZc/JKq8eX3nWiYrLnLt0hVH9PLOipoaE8lsPpA7mV/dPYuvRCzy8ahs1V6xd1a2w1DvjEe0ZHhvGhsxphAU6mLc8l10nKr1eg/KM/JKWK/Ejoz3z41xDQvm1+6Ym8d+zJ5NfUsHDq7ZxudG6J7MLSioIDrAxdmCkJccf0j+U7MwMokIDeHBFHjucoaV6t4KSCiKCHCRFaEgo1SOzJg/ihTktq7o9t72OixYt/1lYWsmkpGgCOnnYz9MGx4SyIXMaMWGBzF+5jfxivR+kt8svrmDykOhrzjbQUxoSqk+4e+JAXpw3leKLzcxf4f3Feuoamth7ssqSS01tJUaHsCFzGnERQSx4ZRt5zmVUVe9zsa6Bg2cuefTvlYaE6jPumpDAk1OC2O9crKey1nuL9ew6UUVjs/GJkABIiAomOzODxOgQHl61nU+PnLe6JNUDO0srMQZSh7pnssj2aEioPmVKnKPVYj15lHtpVberD9FNGeIbIQEQFxnM+iUZDIkJZdHq7Xx8+JzVJaluyi+pwCYweUi0x46hIaH6nBnXxbFiQSpF56qZm5XLeS+s6lZQUkFybJhLiwx5woCIINYtSSd5QDiPrskn56BOkdabFJSUM3ZgJOFB7pirtX0aEqpPumX0AFY9fAOl5bUtq7p5cLEeYwyFpd59iK47+ocHsW5xOqPiwslcW8CW/WesLkl1QWNTMztKKz1+CVNDQvVZXxoZy+pFN3Cq8jKzs3I5XeWZoDh6robymnpSfTQkAPqFBbJucQZjB0bw2GsFvLf3tNUlqWs4cPoStfVNGhJKeVJ6cn9efTSNc5euMDtrK6cq3b+q21bn3UPTRnhmenB3iQoN4NXF6UwYFMUTrxeyeXeZ1SWpTlwd50od5rlBa9CQUIqUoTGsfTSN8up6Zmdt5URFrVv3n1t0gYFRwQyJCXXrfj0hMjiAtY+kMXlwNP+4fgebPmu7NIzyFfklFSREBpMYFezR42hIKAVMHdKP1xanU1XbwOxluZRecE9QGGPIK7pARnJ/xEMPO7lbRHAAax5JI3VoP76TvYM3C09YXZJqR2FJBSnD+nn875WGhFJOkwZH8/riDGrqG5mdtZXi866v6nb0XDXnq+vJSPbsJQF3CwtysGrRDWQk9+d7f/iMjfnHr91IeU1Z1WVOVl72yjiXhoRSrVyfFMW6xRlcaWzmgWVbOXK22qX9bT3qHI9IjnVHeV4VGuhg5cIbuGlkLE+9sYv120qtLkk55Rc7xyM8+BDdVRoSSrUxLrFlsZ5mA3Oytrq0WE9uUTmJUcEMjvH8IkOeEBJoZ/mCVG4dPYBn39zNq1uLrS5J0TJoHRJgZ8zACI8fS0NCqXaMjo8gOzMDmwhzs3I5eLr7QWGMIbeXjUe0JzjATtaCFG4bE8eP3t7Lqk+O9Xhfxhje2nHCY7cb9xUFJRVMHuydySI1JJTqwMi4cLIzM3DYhbnLu7+q25Gz1VyoqScj2bdvfe2KIIedlx5K4c7x8fz7/9vH8r8W9Wg/v/7LYb674TO++dKnbrs5oK+pudLIvrKLpA7zznM3GhJKdSJ5QDjZmdMItNuY181V3a4+H+EPIQEQ6LDx4ryp3H39QH62eT+/zznSrfabPjvFb7Yc5vaxcdTUN/LAsq0UnXNtzKcv+ux4JU3NxmtP8GtIKHUNw2PDyM7MINS5qltXF+vJLbrQq8cj2hNgt/GbOZP5+qREfvGng7z4weEutdtRWsE//+Ez0obF8LsHp5KdmUFDUzOzs3I57MKYT19UUFKBSMtt296gIaFUFwyLDWPD0gyiQwN5aEUeuddYg6G52ZBXVE7GiN49HtEeh93G8w9M4htTBvHcnw/x678c6nT7qssNZL5aQHxkEC/PTyHIYWdMQiQblmYAMHd5rks3B/Q1+SUVjI6LICokwCvH05BQqouS+oWycek0EqKCeXjVNj492vEaDJv3lHGhpp4Z18V5sULvcdhtPHf/JL6VksSv/3KY5947iDHtLw27Lq+Uc5eu8OLcqV+YBXdknOs3B/Q1zc3enyxSQ0KpbkiICmbD0mkMiQll8Zp8Ckr+fvnP+sZmfvGng4xJiOCr1w+0oErvsNuEX3xzInPTBvPih0f4+Z8O/F1QXGlsYtUnx7hpZCyTBkf/3T5GDPjizQEHTnfv5oC+5vDZai7VNXp1skgNCaW6KTY8iNceTScuIoiHX9n+d4PZr+eVUFpey9Mzx2C3+delprZsNuFn917PQxlDWPZRET95Z/8XguLtHac4e+kKmbckd7iP1jcHzMnKZfeJrt8c0NfkO/9T4q07m0BDQqkeiYsM5vUlGUSGBDB/ZR4fOhfruVjXwG8/OMKXRvRn+ugBFlfpHTab8JNZE1h04zBe+eQY/7ZpL8YYmpsNWR8XMXZgJDeP6vyJ8+GxYWxcOo3woJabA9o7Q1NQUFxBbHigVyeL9NxyRkr5uUHRIaxbks6iVdtZtGo7t4weQFxEEOU19Tw7c6zfDVh3RkT48T3jcNiE5R8fo8kYbh0dx5Gz1fx69uQu/V4M6d8y5vPgijzmr9zGyoU3+Pz06t6WX1JBylDPT+rXmp5JKOWCof3D+NN3buGHd49lZ2kFbxSc4OuTErk+Kcrq0rxORPj+V8fy2K0jeC23lG9n7yAxKpi7J3Z9XCYxOoQNSzMYFB3CI6u3k3eNu8j6krOX6igtr/XKfE2taUgo5aJAh43FNyfz0b/M4Id3j+XHXxtndUmWERGevus6npwxktr6JhbfnNztqSPiIoJZtySDxOhgFq3eTn6xXnqClqnBAVK8OB4BLoaEiMSIyPsictj5td3qReQuETkoIkdE5JlrtReR/iLyoYhUi8iLrtSolLf0Cwtk8c3JxIYHWV2KpUSE731lNH/+7i0sunFYj/YxICKI9UsyiI8M5uFV2z9fha0vyy+uINBhY3xipFeP6+qZxDPAFmPMKGCL8/UXiIgd+B0wExgHzBWRcddoXwf8CPhnF+tTSllARBgdH+HStfO4yGDWL8mgf3jLA4wfHDjjxgp7n/ySCiYlRRHksHv1uK6GxCxgjfP7NcC97WyTBhwxxhQZY+qBbGe7DtsbY2qMMX+jJSyUUn1UQlQwbzz2JUbEhbFkbUGfXfyorqGJvaeqSPHyeAS4fndTvDGmDMAYUyYi7T1eOgho/Sd7AkjvRvtOiUgmkAkQHx9PTk5Od3fxuerqapfa9zZ9rb+gfe6tnhhreLFOeOqNXXy25wB3DOt8Sgp/6HNrB8ubaGgyBF06QU7O6Xa38VSfrxkSIvIXIKGdj37QxWO0d77Z/vP7PWCMyQKyAFJTU8306dN7vK+cnBxcad/b9LX+gva5N7ttejNPrCsk+8BZvvnl1E4nuPOXPl+1L+cIcJCFd9/yhalNWvNUn695uckYc7sxZkI7v94GzojIQADn17Pt7OIEMLjV6yTglPP7rrRXSikCHTZ+9cAkBkYF80/rd1B1ucHqkrymoLiC5AFhHQaEJ7k6JrEJWOj8fiHwdjvbbAdGichwEQkE5jjbdbW9UkoBEBkcwAtzp1BWVcf339rd4aSC/sQYQ0FpBSlemhq8LVdD4ufAHSJyGLjD+RoRSRSRzQDGmEbgSeA9YD+w0Rizt7P2zn0UA88DD4vIiVZ3RCml+rCpQ/rxf+4Yzbu7yvhD/gmry/G4o+dqqKxt8Op8Ta25NHBtjLkA3NbO+6eAr7Z6vRnY3NX2zs+GuVKbUsp//cOtI/jwwFn++y+HuG/qIBxeWOvZKlfnsbLizibQJ66VUr2QzSZk3pJMWVUd7+/z7+cn8osriA4NIDk2zJLja0gopXql28bGMyg6hDVbi60uxaOujkfYLJp2XkNCKdUr2W3C/GlDyS0q99vFispr6ik6V+P1+Zpa05BQSvVas1MHE+SwsXZridWleMTVOau8PfNraxoSSqleq19YILMmJ/JW4Um/fG4iv6ScALsw0cKp5zUklFK92oJpw7jc0MQf/HBep4LiCsYnRhEc4N1J/VrTkFBK9WoTBkWRNiyGlz86SkVNvdXluM2VxiZ2nawidah14xGgIaGU8gP/9vXxVNY28B/v7LO6FLfZc/Ii9Y3Nlj1Ed5WGhFKq1xuXGMnjM0by1o6TbNnvH89NXH2IbqqeSSillOuenDGS6+Ij+P5bu6lp6P1zOuUXVzAkJpS4iGBL69CQUEr5hUCHjV/eP5Fzl66wZu8VGpuarS6px4wxFJZWWD4eARoSSik/MjEpmn+5cwzbTjeR+WoBtfWNVpfUIyUXajlfXW/pQ3RXaUgopfzKP0wfwcJxgeQcPMucrFzOXbpidUndlu98iC5FzySUUsr9ZgwJIGt+KofOXOL+lz/lzMU6q0vqloKSciKCHYyOi7C6FA0JpZR/un1cPK8vTufcpSvMXd67zigKSiqYauGkfq1pSCil/FbK0BhWLUqjrLKOectzuVDt+0FRVdvAoTPVPnGpCTQklFJ+Lm14DK88fAPHK2p5cEWez59RFB6/OqmfhoRSSnnFtBH9WbnwBkou1DI7aytlVZetLqlDBcUV2G3CpMHRVpcCaEgopfqIG0fGsvbRNM5evMIDy7ZyvLzW6pLaVVBSwdiBEYQFubS6tNtoSCil+owbhsXw+uJ0Ll5u5IFlWyk+X2N1SV/Q2NTMzuOVpAzxjUtNoCGhlOpjJg2OZv2SDOoampiTletTQbG/7BKXG5pIGWbdIkNtaUgopfqccYmRrFuSQX1TM3OycjnmI0FxdVI/X7mzCTQklFJ91NiBkaxbku4Miq0cPVdtdUkUlFYyMCqYQdEhVpfyOQ0JpVSfNSYhkvVLMmhqNsxelsvB05csraeguNzyqcHb0pBQSvVp1yVEkJ05DZvAnKyt7D1VZUkdpyovc6qqzmeej7hKQ0Ip1eeNjAtn49JphATYmbc8j90nvB8UBT40qV9rGhJKKQUMiw1jw9JpRAQ7eHBFLp8dr/Tq8QtKKggJsDN2YKRXj3stGhJKKeU0OCaU7MwMokIDeGhlHju9GBSFpRVMGhxFgN23fiz7VjVKKWWxpH6hZGdOo19oIPNX5H1+GciTausb2Xvqos9dagINCaWU+juDokPIzswgNiKIBSvzyCu64NHj7SitpKnZcIMPPUR3lUshISIxIvK+iBx2fm03BkXkLhE5KCJHROSZa7UXkTtEpEBEdju/ftmVOpVSqrsSo0PYkJlBQlQwD6/azidHznvsWHnHyrGJ7w1ag+tnEs8AW4wxo4AtztdfICJ24HfATGAcMFdExl2j/Xnga8aY64GFwKsu1qmUUt0WFxlMduY0hsSE8sjq7Xx06JxHjpNXdIHxiVFEBAd4ZP+ucDUkZgFrnN+vAe5tZ5s04IgxpsgYUw9kO9t12N4Ys8MYc8r5/l4gWESCXKxVKaW6bUBEEOszMxgxIJwla/P58OBZt+7/SmMTO45Xkjbc9y41AYgxpueNRSqNMdGtXlcYY/q12eZbwF3GmMXO1/OBdGPMk91o/5gx5vYOasgEMgHi4+NTsrOze9yf6upqwsPDe9y+t+lr/QXtc1/hiT5X1xt+mV/HyUvNPDkliMlx7pnK+1BFE/+ZV8c/TgkiJb7n+3SlzzNmzCgwxqS299k1KxKRvwAJ7Xz0gy4ev71FWruUTCIyHvgv4CsdbWOMyQKyAFJTU8306dO7WNbfy8nJwZX2vU1f6y9on/sKT/X55psaeGhlHr/77CIvzpvAnePb+9HYPXs/PAIcZNE9txATFtjj/Xiqz9e83GSMud0YM6GdX28DZ0RkIIDza3vnYSeAwa1eJwFXLyV12F5EkoC3gAXGmKM96ZxSSrlTVGgAry1OZ3xiFE+8Xsi7u8pc3mfesXKui49wKSA8ydUxiU20DCzj/Pp2O9tsB0aJyHARCQTmONt12F5EooF3gWeNMZ+4WKNSSrlNVEgArz6axpQh0fzj+kL+Z8fJHu+rsamZguJynx2PANdD4ufAHSJyGLjD+RoRSRSRzQDGmEbgSeA9YD+w0Rizt7P2zu1HAj8SkZ3OX3Eu1qqUUm4RERzA6kVppA/vz3c37uSNghM92s++sovU1Df5dEi4NPJijLkA3NbO+6eAr7Z6vRnY3I32PwV+6kptSinlSWFBDlYtuoHFa/L5lzc+wxjD/amDr92wlW3HWhYZ8uWQ0CeulVKqh4ID7KxYmMpNI2N56o+72Jh/vFvt846VM6x/KPGRwR6q0HUaEkop5YLgADvLF7QExdN/3MW6vNIutWtuNmz38fEI0JBQSimXXQ2K6aMH8P23drPi46Jrtvnk6Hkqaxu4cWSsFyrsOQ0JpZRyg+AAO8vmpzJzQgI/fXc/L2w5TGcPKy/7qIj4yCDumuD6sxaepCGhlFJuEuiw8du5U7hvyiCef/8Qv/3gSLvb7TlZxd+OnOeRG4cT5LB7ucrucc9z5UoppQBw2G08d/8kAJ5//xDhQQ4euWn4F7ZZ9tciIoIczE0fYkWJ3aIhoZRSbmazCb/41kRq6hv5j3f2ER7k4IEbWm6PPV5ey7u7TrHk5mQifXDW17Y0JJRSygMcdhsvzJ3C4jX5PPPmLnIOneWbU5P44MBZ7DZh0Y3Dr70TH6AhoZRSHhLksJM1P5Xn/nyQt3acZPPu0wDcn5JEQpTvPhvRmoaEUkp5UEignR/dM46n7xpDzsGzfHToHI/PGGl1WV2mIaGUUl4Q6LDxlfEJfMUN04t7k94Cq5RSqkMaEkoppTqkIaGUUqpDGhJKKaU6pCGhlFKqQxoSSimlOqQhoZRSqkMaEkoppToknc133tuIyDmgxIVdxALn3VROb9DX+gva575C+9w9Q40xA9r7wK9CwlUikm+MSbW6Dm/pa/0F7XNfoX12H73cpJRSqkMaEkoppTqkIfFFWVYX4GV9rb+gfe4rtM9uomMSSimlOqRnEkoppTqkIaGUUqpDGhKAiNwlIgdF5IiIPGN1PZ4gIoNF5EMR2S8ie0Xk2873Y0TkfRE57Pzaz+pa3UlE7CKyQ0Tecb726/4CiEi0iLwhIgecf97T/LnfIvJd59/pPSKyXkSC/a2/IvKKiJwVkT2t3uuwjyLyrPPn2UERudOVY/f5kBARO/A7YCYwDpgrIuOsrcojGoHvGWPGAhnAE85+PgNsMcaMArY4X/uTbwP7W7329/4C/Ab4kzFmDDCJlv77Zb9FZBDwT0CqMWYCYAfm4H/9XQ3c1ea9dvvo/Hc9BxjvbPN758+5HunzIQGkAUeMMUXGmHogG5hlcU1uZ4wpM8YUOr+/RMsPjkG09HWNc7M1wL2WFOgBIpIE3A2saPW23/YXQEQigVuAlQDGmHpjTCX+3W8HECIiDiAUOIWf9dcY81egvM3bHfVxFpBtjLlijDkGHKHl51yPaEi0/KA83ur1Ced7fktEhgFTgDwg3hhTBi1BAsRZWJq7/Rp4Cmhu9Z4/9xcgGTgHrHJeZlshImH4ab+NMSeB54BSoAyoMsb8GT/tbxsd9dGtP9M0JEDaec9v7wsWkXDgj8B3jDEXra7HU0TkHuCsMabA6lq8zAFMBV4yxkwBauj9l1o65LwOPwsYDiQCYSLykLVVWc6tP9M0JFpSdnCr10m0nK76HREJoCUgXjfGvOl8+4yIDHR+PhA4a1V9bnYj8HURKablEuKXReQ1/Le/V50AThhj8pyv36AlNPy137cDx4wx54wxDcCbwJfw3/621lEf3fozTUMCtgOjRGS4iATSMuCzyeKa3E5EhJbr1PuNMc+3+mgTsND5/ULgbW/X5gnGmGeNMUnGmGG0/Jl+YIx5CD/t71XGmNPAcRG5zvnWbcA+/LffpUCGiIQ6/47fRst4m7/2t7WO+rgJmCMiQSIyHBgFbOvpQfSJa0BEvkrL9Ws78Iox5mfWVuR+InIT8DGwm///Gv33aRmX2AgMoeUf3P3GmLYDZL2aiEwH/tkYc4+I9Mf/+zuZlsH6QKAIWETLfwj9st8i8u/AbFru4NsBLAbC8aP+ish6YDot04GfAf4V+B866KOI/AB4hJbfk+8YY/63x8fWkFBKKdURvdyklFKqQxoSSimlOqQhoZRSqkMaEkoppTqkIaGUUqpDGhJKKaU6pCGhlFKqQ/8fef9Abc4Fc/cAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.plot(y.T[100])\n", "plt.grid()" ] }, { "cell_type": "code", "execution_count": 86, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([0.00000000e+00, 3.14158749e-03, 6.28314397e-03, ...,\n", " 6.28314397e-03, 3.14158749e-03, 1.22464680e-16])" ] }, "execution_count": 86, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y0" ] }, { "cell_type": "code", "execution_count": 87, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 87, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAqbUlEQVR4nO3dd3xUZb7H8c8vvZCQBBJKCoEUeg9FFKSKKIpdsK+FRUXUXd1l71p2716vuuLqYkNEXSssVtCFRUEQGyVIEYSQEAgJkEJNI/25f2T0ZjGQSZjJmfJ7v155kTlzJvM9gXx5cspzxBiDUkop9+djdQCllFKOoYWulFIeQgtdKaU8hBa6Ukp5CC10pZTyEH5WvXH79u1NYmKiVW+vlFJuadOmTYeNMdGNPWdZoScmJpKenm7V2yullFsSkZzTPae7XJRSykNooSullIfQQldKKQ+hha6UUh5CC10ppTxEk4UuIq+JSKGIbD/N8yIic0UkS0S2icggx8dUSinVFHtG6P8ALjzD85OAFNvHdOCls4+llFKquZo8D90Ys1ZEEs+wyhTgTVM/D+86EYkQkU7GmEOOCqmUo9XVGfKLK8gvruBIaRWHSyspqaimps5QW2uoMxAa6EubQD/aBPkR3SaQuKgQOoYH4esjVsdXqlGOuLAoFsht8DjPtuwXhS4i06kfxZOQkOCAt1aqadW1dew8VMymnGNszT1OZmEp2UVlnKyubfbX8vMR4qNC6NU5nN6dw+kb25ZBCZGEBlp2jZ5SP3PEv8LGhiuN3jXDGDMfmA+Qlpamd9ZQTpN7tJwvdhWyalchG/YeoaK6DoCO4UGkdgxjWNd2dIsOJTYymPahgbRrE0B4sD9+PoKfjyAilFfVUFpZQ0lFDQXFFeQdO0nesXL2FJaxLe84/9pWP2bx8xEGJURybnJ7xvaIoU9sOCI6iletzxGFngfEN3gcBxx0wNdVqlnyT1Tw8ZYDfLz5ALvySwDoFh3K1CEJpCVGMighks4RwXZ/vbAgf8KC/OnUFlI7hP3i+ePlVWzLO8G3e47wTdZhnl21m2dW7iYhKoSL+3Xikn6d6dU53GHbp1RTHFHoS4GZIrIIGAac0P3nqrXU1hk+/7GAd9bn8HXWYYyBwV0ieejinozr2YGu7UOd9t4RIQGMSo1mVGr9PElHy6pY+WMBn/5wiPlrs3lpzR76x7XlumEJXNK/MyEBultGOZc0dU9REVkIjAbaAwXAo4A/gDFmntT/bvk89WfClAO/MsY0OetWWlqa0cm5VEsVV1SzeGMu//h2H3nHThIbEcyVg2K5fFCcU0vcXsfKqli69SDvrM9hd0EpYYF+XDc8gdvO60pMWJDV8ZQbE5FNxpi0Rp+z6ibRWuiqJUoqqvnHN/t45atsiitqGJoYxa3nJTK+Zwf8fF3vOjljDOk5x3jj230s++EQfr4+XD04jhnnJxEfFWJ1POWGzlTo+jugcgsnq2p5/du9zF+bzfHyasb37MCsccn0i4uwOtoZiQhDEqMYkhjFvsNlvLw2m/fS81icnsuNwxO5Z2wykaEBVsdUHkJH6MqlGWP4dNshHl+2k4MnKhjTPZr7xqfSPz7C6mgtln+igmdX7mZxei6hAX7cOSaJW8/tSpC/r9XRlBvQXS7KLe04eII/f/IjG/YepWencP50SS+GdWtndSyHySwo4cl/72LlzkIS24Xwl8v6MDKl0RvRKPUzLXTlViqqa3l2ZSavfJVNeJAfD0zsztQhCR57heY3WYd56OPt7D1cxqX9O/PQ5J564FSdlha6chubco7y4PvbyC4q45q0OP54US/ahvhbHcvpKqprmfflHl5cvYcgfx/+clkfLu3fWS9QUr+gha5cXmVNLXNWZLDg6710bhvME1f29crdD9lFpTzw3la+33+ci/t14n+m9NGDpuo/6FkuyqVlF5Uya9Fmth8o5vphCfzXRT29dm6UbtFtWPzrc3h5bTbPfL6bjXuPMufq/j9fvKTUmbjeibvKaxhj+GBTHpOf+5q8YyeZf+NgHru8r9eW+U/8fH24e0wyH999Lm2D/bn59Q387fPd1Nbp9EfqzLTQlSUqa2r5w4c/8Nv3ttI3ti3L7x3JBb07Wh3LpfSJbcvSmedxxcA45q7K5ObXNnC4tNLqWMqFaaGrVldQXMHU+etYtDGXu0Yn8e4dw+nU1v5Js7xJcIAvc67ux5NX9mXjvqNcPPcrNu8/ZnUs5aK00FWr2pRzlMnPfU1GfgkvXj+I313Yw2NPR3QUEeHaIQl8dNe5BPj5cO38dSzZcsDqWMoFaaGrVvPR5jymzl9HSIAvH911Lhf17WR1JLfSq3M4S+4+jwFxEdy7aAt/+3w3dbpfXTWgha6czhjD819kcv8/t5LWJYqld59H946/nF9cNS0qNIC3bx/G1YPr96vfs3AzFS2485LyTN59OoFyupraOh76eDuLNuZy+cBYnryyHwF+Oo44GwF+Pvz1qn6kdgjjf5fvpLCkggU3D6FtsOdfgKXOTH+ylNOUV9Vw+5vpLNqYy8wxyfztmv5a5g4iItwxqhvPTRvIltzjXPvydxQUV1gdS1lMf7qUUxRXVHPTqxtYu7uI/728Lw9M7K6XsTvB5H6def2WoeQeLefKl74lu6jU6kjKQlroyuGOlVVxw4L1bMk9znPTBnHdsASrI3m081Las3D6cE5W1XL1vO/YcfCE1ZGURbTQlUMVltSfY74rv4T5Nw3m4n56Jktr6BcXwXszziHQz4frF6xn+wEtdW+kha4c5tCJk0x9eR25x8r5xy1DGNujg9WRvEq36DYsmn4OoQF+XPfKOrblHbc6kmplWujKIQqLK5g2fx1FJZW8ddtQRiS3tzqSV0poF8Ki6cMJD/bn+gXr9apSL6OFrs7a4dJKrluwnqKSSv5x61AGd4myOpJXi48K4Z+/PofIkABuenUDW3OPWx1JtRItdHVWjpfXHwDNO1bOa7cMYXCXSKsjKSA2IphF04fTNqR+tsaM/BKrI6lWoIWuWqy4opobX91A9uEyFtw0xKPu9+kJOkcE887twwjw9eGGV9eTc6TM6kjKybTQVYtUVNdy6+sb2ZVfzLwbBnFeiu4zd0Vd2oXy9u3DqKmt4/oF6zl04qTVkZQTaaGrZquprWPmu5vZtP8Yz147UM9mcXGpHcJ449ahHC+v5oYF6zmic6p7LC101SzGGB5esp2VOwv486W99TxzN9EvLoLXbhlC3rGT3PpGOuVVNVZHUk6gha6a5ZmVmSzcUD83y03nJFodRzXD0K5RzJ02kG15x5m1cIve0s4DaaEru729Loe5qzK5Ni2e316QanUc1QITe3fkz5f2ZuXOAh5duh1jtNQ9iU6fq+yyelchjyzZzrgeMTx2eR+daMuN3XROIgeOneTltdnERoRw5+gkqyMpB9FCV03alV/MPQs307NTOM9dNxA/X/3Fzt39/sIeHDxRwZP/3kXniCCmDIi1OpJyAC10dUaHSyu57R/phAb68urNQwgJ0H8ynsDHR5hzdT8Kiit48P1tJESFMDBBLwpzd3YNtUTkQhHJEJEsEZndyPNtReQTEdkqIjtE5FeOj6paW0V1LdPfTOdIWSULbhpCx7ZBVkdSDhTo58u8GwbTITyQX7+1ifwTeoMMd9dkoYuIL/ACMAnoBUwTkV6nrHY38KMxpj8wGnhaRAIcnFW1ImMMv/9gG9/vP86z1w6gb1xbqyMpJ4gKDWDBTUMoq6xh+lvpnKzS+5O6M3tG6EOBLGNMtjGmClgETDllHQOESf2RsjbAUUBPdHVjL6zOYsmWgzw4sTsX9tFzzT1Z945hPDt1ID8cOMHvPtimZ764MXsKPRbIbfA4z7asoeeBnsBB4AfgXmNM3alfSESmi0i6iKQXFRW1MLJyttUZhTz9+W4uG9CZu/QMCK8woVcHHpzYnU+2HuSF1VlWx1EtZE+hN3Z+2qn/hU8EtgCdgQHA8yIS/osXGTPfGJNmjEmLjo5uZlTVGnKOlHHvws307BjO41f009MTvcid5ydx2YDOzPlsNyt/LLA6jmoBewo9D4hv8DiO+pF4Q78CPjT1soC9QA/HRFStpbyqhl+/tQkfH+HlGwcTHOBrdSTVikSEJ67sR5/YcO5fvIX9R8qtjqSayZ5C3wikiEhX24HOqcDSU9bZD4wDEJEOQHcg25FBlXPVHwT9gYyCEuZOHUh8VIjVkZQFgvx9een6wfiIMOPtTVRU60FSd9JkoRtjaoCZwApgJ7DYGLNDRGaIyAzban8BRojID8Aq4PfGmMPOCq0c79Wv9/LJ1oM8cEF3RqXq7jBvFh8VwrPXDmBnfjEPfazTA7gTu64SMcYsA5adsmxeg88PAhc4NppqLRv3HeXx5bu4sHdHPQiqABjTI4Z7xqYwd1UmgxIiuW5YgtWRlB30Gm4vd7Ssinve3Ux8ZDBPXa0HQdX/u3dcCiNT2vOnpTvYlnfc6jjKDlroXqyuzvDbxVs4WlbF89cNIizI3+pIyoX4+ghzpw4kOiyQO9/+nhMnq62OpJqghe7FXvkqm9UZRTw8uSd9YvVKUPVLkaEBPH/dQAqKK/jDh3rRkavTQvdSm3KO8dcVGVzUtyM3DO9idRzlwgYmRPLAxO4s+yGfhRtym36BsowWuhc6Xl7FrIWbiY0I5okrdb+5atr0kd0YmdKeP3+yg4z8EqvjqNPQQvcyxhgeeG8bhSUVPH/dQMJ1v7myg4+P8LdrBhAW5M/Md7/XSbxclBa6l3l7/X5W7ixg9qSe9IuLsDqOciPRYYE8c21/MgtL+e9Pd1gdRzVCC92LZBWW8ti/fmRUajS3nptodRzlhkamRHPn6CQWbsjlk62nzgCirKaF7iWqauq475+bCfb3Zc5Vut9ctdxvJqQyMCGCP370A4dOnLQ6jmpAC91LPLNyN9sPFPPElf2ICdc7D6mW8/f14ZlrBlBda3jgva3U1empjK5CC90LrMs+wrwv9zB1SDwTe3e0Oo7yAIntQ3l4ci++yTrCG9/tszqOstFC93AnTlbzm39uoUtUCA9PPvXOgUq13LSh8YztEcMTy3eRVainMroCLXQP98iS7RSUVPLMtQMIDbRrLjal7FI/f3pfQgJ8ue+fW6iq+cVNylQr00L3YJ9uO8iSLQeZNTaFgQmRVsdRHigmLIjHr+jL9gPFPPdFptVxvJ4Wuoc6XFrJI0t20C+uLXeP0SlxlfNc2KcTVw2O44XVWWzKOWZ1HK+mhe6BjDE8/PF2SitqmHN1f/x89a9ZOdejl/SiU9tgfrt4i15FaiH9SfdAn247xPLt+dw3IYXUDmFWx1FeICzInzlX92ffkXLmfJZhdRyvpYXuYYpKKnlkyXb6x0cwfWQ3q+MoL3JOUjtuGJ7Aa9/s1V0vFtFC9yA/7Wopq6xlzlX9dFeLanWzJ/Wkc9tgfvf+Vr3BtAX0J96DfLLtEP/ekc/9E1JJ0V0tygJtAv14/Iq+7Ckq49mVetZLa9NC9xCFJRU/72q5Y2RXq+MoLzYqNZpr0+KZv3YPW3OPWx3Hq2ihe4hHl+ygvKqWp6/WXS3Ken+c3JOYsCAefH8rlTW666W16E++B1ixI5/l2/O5d1wKyTG6q0VZLzzIn8ev6MvuglKe/yLL6jheQwvdzZVUVPPIku306BjG9FF6VotyHWN6xHDFoFheXLOHnYeKrY7jFbTQ3dxf/51BYUklT1zZD3/d1aJczMMX9yIi2J8/fPgDtTrNrtNpA7ixTTlHeXt9DreMSGRAfITVcZT6hcjQAB6e3Istucd5e12O1XE8nha6m6qsqWX2Bz/QKTyI317Q3eo4Sp3WlAGdGZnSnqdWZJB/osLqOB5NC91NzVuTTWZhKf9zeR/a6LS4yoWJCP9zWR+qa+t4dOl2q+N4NC10N5RVWMILq7OY3K8TY3t0sDqOUk3q0i6U+8ansmJHASt25Fsdx2NpobuZujrDf324neAAXx69pLfVcZSy2+0ju9KjYxiPLtlBSUW11XE8kl2FLiIXikiGiGSJyOzTrDNaRLaIyA4R+dKxMdVP3t+Ux4Z9R/mvi3oQHRZodRyl7Obv68PjV/SloKSCpz/bbXUcj9RkoYuIL/ACMAnoBUwTkV6nrBMBvAhcaozpDVzt+KjqWFkVjy/fSVqXSK4eHG91HKWabWBCJDcO78Ib3+1ji04L4HD2jNCHAlnGmGxjTBWwCJhyyjrXAR8aY/YDGGMKHRtTAfx1RQbFFTX85bI++PiI1XGUapEHJ3YnJiyQP36k56Y7mj2FHgvkNnicZ1vWUCoQKSJrRGSTiNzU2BcSkekiki4i6UVFRS1L7KU27z/Goo37uWVEIj07hVsdR6kWCwvy56GLe7HjYDHvrtdz0x3JnkJvbCh46n+rfsBg4GJgIvCwiKT+4kXGzDfGpBlj0qKjo5sd1lvV1hkeXrKd6DaB3Dc+xeo4Sp21yf06MSKpHU+tyOBwaaXVcTyGPYWeBzTcYRsHHGxknX8bY8qMMYeBtUB/x0RU76zPYfuBYh6e3IuwIH+r4yh11kSE/57Sm5PVtTy5fJfVcTyGPYW+EUgRka4iEgBMBZaess4SYKSI+IlICDAM2OnYqN6pqKSSp1ZkcF5yeyb362R1HKUcJjkmjNvO68Z7m/LYlHPU6jgeoclCN8bUADOBFdSX9GJjzA4RmSEiM2zr7AT+DWwDNgALjDF6SZgDPL5sJxXVtfx5Sm9E9ECo8iz3jE2mU9sgHvp4BzW1dVbHcXt2nYdujFlmjEk1xiQZYx6zLZtnjJnXYJ2njDG9jDF9jDHPOimvV1mXfYQPNx9g+qhuJEW3sTqOUg4XGujHw5N7sfNQsU7e5QB6paiLqq6t45El24mNCGbmGD0QqjzXpD4dGZnSnqc/201RiR4gPRta6C7qjW/3sbuglD9d2pvgAF+r4yjlNCLCny7tTUVNLY8v10NvZ0ML3QUdLq3k76syGZUazfieMVbHUcrpkqLbMH1UNz78/gAb9+kB0pbSQndBT3+WwcmqWh6Z3EsPhCqvcfeY+gOkf1q6Q68gbSEtdBez/cAJFm3M5eYRiSTH6IFQ5T1CAvyYPakHOw4W8/6m3KZfoH5BC92FGGP409IdRIUEMGucHghV3ufS/p1J6xLJUysyKNYpdptNC92FfLLtEOk5x3hgYnfaBusVocr7iAiPXtKbI2VVPP9FltVx3I4Wuosor6rh8WU76d05nGvSdGpc5b36xrXl6sFxvP7NXrKLSq2O41a00F3EvDV7OHSigj9d2htfnRpXebkHJnYn0M+Xx/6lpzE2hxa6C8g9Ws7La7O5pH9nhiRGWR1HKcvFhAUxa1wyq3YVsiZDb69gLy10F/D48p2IwB8m9bA6ilIu45YRXenaPpS/fPoj1TrPi1200C323Z4jLPshn7tGJ9M5ItjqOEq5jAA/Hx66uCd7isp48zud58UeWugWqq0z/PenPxIbEcz0Ud2sjqOUyxnbI4ZRqdE8u3I3R/RGGE3SQrfQB9/nsfNQMb+f1IMgf52vRalTiQiPTO5JeVUtz67MtDqOy9NCt0h5VQ1zVmQwID6CS/TGFUqdVnJMGNcNTeDdDfvJKiyxOo5L00K3yPy12RSWVPLQxT11vhalmnDf+BSC/X15Qm9Xd0Za6BYoKK7g5S+zmdSnI2l6mqJSTWrXJpC7xiSxcmch3+45bHUcl6WFboGnP8ugpq6O2XqaolJ2u/XcrsRGBPPYv3ZSp7MxNkoLvZX9eLCY9zblcdM5iXRpF2p1HKXcRpC/Lw9O7M6Og8V8tPmA1XFckhZ6KzLG8L/LdhIe5M89Y5OtjqOU27m0f2f6xbVlju2eAeo/aaG3ojUZRXyddZhZ41KICAmwOo5SbsfHR/jjRT05dKKCV7/OtjqOy9FCbyU1tXU8tmwnie1CuHF4F6vjKOW2hnVrxwW9OvDSmj16U+lTaKG3kkUbc8kqLGX2pB4E+Om3XamzMXtSDypr6nhm5W6ro7gUbZZWUFJRzTOf72ZoYhQTe3e0Oo5Sbq9bdBtuGN6FRRv2k1mgFxv9RAu9FbyyNpsjZVX8l15EpJTDzBqXQmigH4/rxUY/00J3ssKSCl75ai8X9+vEgPgIq+Mo5TGiQgOYOSaZL3bpxUY/0UJ3srmrMqmureOBC7pbHUUpj3PziEQ6tw3iyeW7MEYvNtJCd6K9h8tYuCGXaUMT6NpeLyJSytGC/H25b0IqW/NOsHx7vtVxLKeF7kRzPssg0M+He8bpRURKOcuVg+JI7dCGOSsyvP7ORlroTrI19zj/2naI28/rSkxYkNVxlPJYvj7CgxN7kH24jMXpuVbHsZQWuhMYY3hi+S6iQgO4Q+9EpJTTje8ZQ1qXSP6+MpPyqhqr41jGrkIXkQtFJENEskRk9hnWGyIitSJyleMiup+1mYf5LvsI94xNJizI3+o4Snk8EWH2pB4UllTy+jf7rI5jmSYLXUR8gReASUAvYJqI9DrNek8CKxwd0p3U1dWPzuOjgrluWILVcZTyGmmJUYzv2YF5a/ZwrKzK6jiWsGeEPhTIMsZkG2OqgEXAlEbWuwf4ACh0YD6388m2g+w8VMxvJ3Qn0E/vE6pUa/rdhd0pq6rhxTVZVkexhD2FHgs0PNKQZ1v2MxGJBS4H5p3pC4nIdBFJF5H0oqKi5mZ1eZU1tTy1IoOencK5tH9nq+Mo5XVSO4Rx5aA43vg2hwPHT1odp9XZU+iNXat+6hn8zwK/N8accYJiY8x8Y0yaMSYtOjrazoju4931+8k7dpLZk3rg46OX+CtlhfsnpILAM59738Rd9hR6HhDf4HEccPCUddKARSKyD7gKeFFELnNEQHdRUlHNc19kMSKpHaNS2lsdRymv1TkimFtGJPLB93lk5HvXxF32FPpGIEVEuopIADAVWNpwBWNMV2NMojEmEXgfuMsY87Gjw7qyV9Zmc7Ssit9f2EMn4FLKYneNTqJNoB9PrfCuibuaLHRjTA0wk/qzV3YCi40xO0RkhojMcHZAd3C4tJIFX+/lor4d6a8TcClluYiQAO4cncTKnYVs2HvU6jitxq7z0I0xy4wxqcaYJGPMY7Zl84wxvzgIaoy5xRjzvqODurIXV++horqW30zQCbiUchW/GtGVmLBA5qzI8JqJu/RK0bN06MRJ3l6fw5WD4kiOaWN1HKWUTXCAL/eMTWbDvqOszfSO6XW10M/S3FVZGGOYNS7F6ihKqVNcOySBuMhgnv7MO0bpWuhnYZ9tMqBpQxOIjwqxOo5S6hQBfj7cOy6FbXknWLGjwOo4TqeFfhaeXbkbf19h5hidHlcpV3X5wFi6RYfy9GcZ1NZ59ihdC72FMvJLWLL1IDePSCQmXKfHVcpV+fn68JsJqWQWlrJ06wGr4ziVFnoL/e3zDNoE+DFjVJLVUZRSTbioTyd6dgrnmc8zPfomGFroLbA19zgrdhRw28iuRIYGWB1HKdUEHx/hwYmp7D9aznvpeVbHcRot9BaY81kGkSH+3HZeV6ujKKXsNKZ7DIMSIpi7KpOK6jNOO+W2tNCbaX32Eb7KPMydo5P05hVKuRER4YGJ3ckvruDtdTlWx3EKLfRmMMYw57MMYsICuemcRKvjKKWaaURSe85NbsdLa/ZQVul5t6rTQm+GNbuL2LjvGPeMTSbIX29eoZQ7euCC7hwpq+L1b/ZaHcXhtNDtZIzh6c8yiIsM5tohems5pdzVwIRIxveM4eW12Zwor7Y6jkNpodvp39vz2X6gmPvGpxLgp982pdzZbyZ0p6Sihvlf7bE6ikNpM9mhts7w9Oe7SYoO5fKBsU2/QCnl0np1Dmdyv068/s0+ikoqrY7jMFrodliy5QBZhaX8ZkJ3fPXWckp5hPsnpFJRXetRN5TWQm9CTW0dc1dl0rNTOJP6dLQ6jlLKQZKi23DFoDjeWb+fguIKq+M4hBZ6Ez7afIB9R8q5f3yK3vhZKQ8za2wKtXWGl9Z4xr50LfQzqK6tY+4XmfSJDWdCrw5Wx1FKOVhCuxCuHhzHu+v3c+jESavjnDUt9DP48Ps8co+e5P7xqXrjZ6U81N1jkqkzhhdXu/8oXQv9NKpq6pi7Kov+cW0Z2yPG6jhKKSeJjwrhmiHxLNq4nwPH3XuUroV+Gu9vyuPA8ZPcp6NzpTze3bab1Lyw2r3PeNFCb0RlTS0vrM5iQHwEo7tHWx1HKeVksRHBTB2SwOKNueQeLbc6TotpoTdicXr96Pz+CTo6V8pb3DUmCR8Rnv/CfUfpWuinqKiu5cXVWQzuEsmolPZWx1FKtZJObYO5blgC73+fR86RMqvjtIgW+in+uTGXQycq9MwWpbzQnaOT8PMRnnPTUboWegM/XQY8NDGKc5PbWR1HKdXKOoQHcf2wLny0+QB7D7vfKF0LvYF31++noLiS+yak6OhcKS81Y3Q3/H2F51ZlWh2l2bTQbSqqa3npyz0M6xrFiCTdd66Ut4oJC+LG4V34eMsB9hSVWh2nWbTQbd5el0NRSSX3T0i1OopSymK/Pj+JQD9f5rrZKF0LHSivqmHel3sYkdSO4d1037lS3q59m0BuHpHI0q0HySossTqO3ewqdBG5UEQyRCRLRGY38vz1IrLN9vGtiPR3fFTneXtdDodLq3R0rpT62fRR3Qjx9+XZle4zSm+y0EXEF3gBmAT0AqaJSK9TVtsLnG+M6Qf8BZjv6KDOUlZZw7wvsxmZ0p4hiVFWx1FKuYio0ABuOTeRf/1wiIx89xil2zNCHwpkGWOyjTFVwCJgSsMVjDHfGmOO2R6uA+IcG9N53vwuh6NlVdw3XkfnSqn/dMfIboQG+PH3VbutjmIXewo9Fsht8DjPtux0bgOWN/aEiEwXkXQRSS8qKrI/pZOUVtbw8to9nJ8azeAukVbHUUq5mIiQAG49N5FlP+SzK7/Y6jhNsqfQGzsh2zS6osgY6gv99409b4yZb4xJM8akRUdbP+nVG9/u43h5te47V0qd1q3ndaVNoB/PrXL9q0ftKfQ8IL7B4zjg4KkriUg/YAEwxRhzxDHxnKessoYFX2Uzuns0A+IjrI6jlHJRESEB3DIikWXbD7G7wLX3pdtT6BuBFBHpKiIBwFRgacMVRCQB+BC40RjjFjub3lqXw7HyamaNS7E6ilLKxd12XldC/F3/vPQmC90YUwPMBFYAO4HFxpgdIjJDRGbYVnsEaAe8KCJbRCTdaYkdoLyqhlfW1p/ZMihB950rpc4sMjSAm0bUn/Hiyuel23UeujFmmTEm1RiTZIx5zLZsnjFmnu3z240xkcaYAbaPNGeGPlvvrt/PkbIq7tXRuVLKTneM7Eawv69Lz8TodVeKnqyqZd6X2Zyb3I40Pe9cKWWnqNAAbjynC59sPeiyc7x4XaEv3LCfw6WVzBqro3OlVPPcMbIbgX6+LntXI68q9IrqWubZZlQcpnO2KKWaqX2bQG4YnsCSLa45X7pXFfri9FwKSyp137lSqsWmj0oiwM/HJUfpXlPolTW1vLRmD0MSIzknSUfnSqmWiQ4L5Pph9fOlu9q9R72m0N9Lz+PQiQpmjdO7ESmlzs6vR3XDz0dcbpTuFYVeVVPHS2v2MDAhgvOS9W5ESqmzExMexLShCXy4+QC5R8utjvMzryj0D7/P48Dxk9yro3OllIPcOToJXx/hhdWuM0r3+EKvrq3j+dVZ9I9ry/mp1k8IppTyDB3Cg5g2JJ73N+W5zCjd4wv9o80HyDt2UvedK6UcbsboJHxEeOnLPVZHATy80Gtq63hhdRZ9YsMZ2yPG6jhKKQ/TqW0w1wyJ4730XA4cP2l1HM8u9KVbD5JzpJxZY3V0rpRyjjtHJwPw0hrr96V7bKHX1hme/yKLnp3CmdCrg9VxlFIeKjYimKvT4lm8MY9DJ6wdpXtsoX+67SDZh8uYNTZZR+dKKae6a3QSdcYwb421+9I9stBr6wzPfZFF9w5hTOzd0eo4SikPFxcZwlWD41i4MZeC4grLcnhkoS/ffoiswlLuGZeMj4+OzpVSznf3mGTq6gwvWThK97hCr6szzF2VSXJMGyb16WR1HKWUl4iPCuGKQbEs3LCfQotG6R5X6Ct25LO7oJR7xibjq6NzpVQruntMMjV1hpfXZlvy/h5V6HV1hr+vyqRb+1Am9+tsdRyllJfp0i6UywbE8s76HApLWn+U7lGFvnJnAbvyS5ipo3OllEVmjk2mqqaOBV/tbfX39phCN6Z+dN6lXQiX9tfRuVLKGl3bh3Jp/8689V0OR0orW/W9PabQv9hVyI6Dxdw9Jhk/X4/ZLKWUG5o5NpmKmlpe/bp1R+ke0XzG1J/ZEh8VzOUDY62Oo5TycskxYVzctxNvfLuP4+VVrfa+HlHoX+4uYmveCe4enYy/js6VUi5g5thkyqpqee2bfa32nm7ffj/tO4+NCOaKQXFWx1FKKQB6dAznwt4def2bvZw4Wd0q7+n2hf511mE27z/OXWPq78StlFKuYubYZEoqanjj232t8n5u3YDGGP6+MpNObYO4arCOzpVSrqVPbFvG94zh1a/3UlpZ4/T3c+tC/y77COk5x7hzdBKBfr5Wx1FKqV+4Z2wKJ05W8+Z3+5z+Xm5d6H9fmUmH8ECuSYu3OopSSjWqf3wE56dGs+CrvZRXOXeU7raFvi77COv3HmXG+UkE+evoXCnlumaNS+FoWRXvrNvv1Pdx20J/7otM2rcJZNrQBKujKKXUGQ3uEsm5ye14eW02FdW1TnsfuwpdRC4UkQwRyRKR2Y08LyIy1/b8NhEZ5Pio/y9931G+yTrCjPO76ehcKeUWZo1N4XBpJQs3OG+U3mShi4gv8AIwCegFTBORXqesNglIsX1MB15ycM7/8PdVmbRvE8D1w7o4822UUsphhnVrx9CuUcz7co/TRun2jNCHAlnGmGxjTBWwCJhyyjpTgDdNvXVAhIg45e4S3+8/xleZh7ljZDeCA3R0rpRyH/eOS6GguJL3NuU55evbU+ixQG6Dx3m2Zc1dBxGZLiLpIpJeVFTU3KwAGAMjU9pzw3AdnSul3MuIpHZc0r8zEcH+Tvn6fnas09jE4qYF62CMmQ/MB0hLS/vF8/YY3CWSt24b1pKXKqWUpUSE56YNdNrXt2eEngc0PNE7DjjYgnWUUko5kT2FvhFIEZGuIhIATAWWnrLOUuAm29kuw4ETxphDDs6qlFLqDJrc5WKMqRGRmcAKwBd4zRizQ0Rm2J6fBywDLgKygHLgV86LrJRSqjH27EPHGLOM+tJuuGxeg88NcLdjoymllGoOt71SVCml1H/SQldKKQ+hha6UUh5CC10ppTyE1B/PtOCNRYqAnBa+vD1w2IFx3IFus3fQbfYOZ7PNXYwx0Y09YVmhnw0RSTfGpFmdozXpNnsH3Wbv4Kxt1l0uSinlIbTQlVLKQ7hroc+3OoAFdJu9g26zd3DKNrvlPnSllFK/5K4jdKWUUqfQQldKKQ/hdoXe1A2r3ZGIxIvIahHZKSI7RORe2/IoEflcRDJtf0Y2eM0fbN+DDBGZaF36syMiviKyWUQ+tT326G0WkQgReV9Edtn+vs/xgm2+3/bveruILBSRIE/bZhF5TUQKRWR7g2XN3kYRGSwiP9iemysijd086PSMMW7zQf30vXuAbkAAsBXoZXUuB2xXJ2CQ7fMwYDf1N+T+KzDbtnw28KTt8162bQ8Eutq+J75Wb0cLt/03wLvAp7bHHr3NwBvA7bbPA4AIT95m6m9FuRcItj1eDNziadsMjAIGAdsbLGv2NgIbgHOovwvccmBSc3K42wjdnhtWux1jzCFjzPe2z0uAndT/IEyhvgCw/XmZ7fMpwCJjTKUxZi/189APbdXQDiAiccDFwIIGiz12m0UknPof/FcBjDFVxpjjePA22/gBwSLiB4RQfzczj9pmY8xa4Ogpi5u1jSLSCQg3xnxn6tv9zQavsYu7FbpdN6N2ZyKSCAwE1gMdjO3OT7Y/Y2yrecr34Vngd0Bdg2WevM3dgCLgddtupgUiEooHb7Mx5gAwB9gPHKL+bmaf4cHb3EBztzHW9vmpy+3mboVu182o3ZWItAE+AO4zxhSfadVGlrnV90FEJgOFxphN9r6kkWVutc3Uj1QHAS8ZYwYCZdT/Kn46br/Ntv3GU6jftdAZCBWRG870kkaWudU22+F023jW2+5uhe6xN6MWEX/qy/wdY8yHtsUFtl/DsP1ZaFvuCd+Hc4FLRWQf9bvOxorI23j2NucBecaY9bbH71Nf8J68zeOBvcaYImNMNfAhMALP3uafNHcb82yfn7rcbu5W6PbcsNrt2I5kvwrsNMb8rcFTS4GbbZ/fDCxpsHyqiASKSFcghfqDKW7DGPMHY0ycMSaR+r/HL4wxN+DZ25wP5IpId9uiccCPePA2U7+rZbiIhNj+nY+j/hiRJ2/zT5q1jbbdMiUiMtz2vbqpwWvsY/XR4RYcTb6I+rNA9gB/tDqPg7bpPOp/tdoGbLF9XAS0A1YBmbY/oxq85o+270EGzTwS7mofwGj+/ywXj95mYACQbvu7/hiI9IJt/jOwC9gOvEX92R0etc3AQuqPEVRTP9K+rSXbCKTZvk97gOexXc1v74de+q+UUh7C3Xa5KKWUOg0tdKWU8hBa6Eop5SG00JVSykNooSullIfQQldKKQ+hha6UUh7i/wDgsYRFCBs+OwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.plot(y0)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "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": 4 }