Skip to content

Instantly share code, notes, and snippets.

@jrgriffiniii
Created April 26, 2020 04:58
Show Gist options
  • Save jrgriffiniii/e477f323b13c97844d01f2792e969c21 to your computer and use it in GitHub Desktop.
Save jrgriffiniii/e477f323b13c97844d01f2792e969c21 to your computer and use it in GitHub Desktop.
Example Jupyter Notebook demonstrating a simple ODE
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import scipy.integrate as spi\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"m = 1. # particle's mass\n",
"k = 1. # drag coefficient\n",
"g = 9.81 # gravity acceleration"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"# The initial position is (0, 0).\n",
"v0 = np.zeros(4)\n",
"# The initial speed vector is oriented\n",
"# to the top right.\n",
"v0[2] = 4.\n",
"v0[3] = 10."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"def func1(v, t0, k):\n",
" # v has four components: v=[u, u'].\n",
" u, udot = v[:2], v[2:]\n",
" # We compute the second derivative u'' of u.\n",
" udotdot = -k / m * udot\n",
" udotdot[1] -= g\n",
" # We return v'=[u', u''].\n",
" return np.r_[udot, udotdot]"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(0.0, 6.0)"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAD8CAYAAACvt3fBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdd3xUVfr48c+dng4pJIQQQBJCC4RiCCxdkaKg66oUBdGlSVNRV1fX/e5PXVfXgitNKaIogmVdUaSIShEEQpEaGAg9kEASAimTZObOvb8/xiAhk5Ayk5kk5/168YrMubn3IcI8c+4953kkVVURBEEQBKF2aTwdgCAIgiA0RCIBC4IgCIIHiAQsCIIgCB4gErAgCIIgeIBIwIIgCILgASIBC4IgCIIHVCoBS5LUSJKkLyVJOipJ0hFJknq6OzBBEARBqM90lTzuP8A6VVXvkyTJAPi6MSZBEARBqPekmxXikCQpCNgH3KKKqh2CIAiC4BKVmQG3AjKBpZIkdQb2AI+rqlpw/UGSJE0CJgH4+fl1a9u2ratjFQRBEASvtGfPnixVVcOq8j2VmQF3B3YAf1BVdackSf8BclVVfbG87+nevbu6e/fuqsQhCIIgCHWWJEl7VFXtXpXvqcwirDQgTVXVnb/9/kuga1WDEwRBEAThdzdNwKqqZgDnJEmK++2l24AUt0YlCIIgCPVcZVdBzwCW/7YC+iTwiPtCEgRBEIT6r1IJWFXVfUCV7m0LgiAI9YvNZiMtLY2ioiJPh+IxJpOJqKgo9Hp9jc9V2RmwIAiC0MClpaUREBBAy5YtkSTJ0+HUOlVVyc7OJi0tjVatWtX4fKIUpSAIglApRUVFhISENMjkCyBJEiEhIS67AyASsCAIglBpDTX5lnDln18kYEEQBEHwAJGABUEQBLeQ7QoWq4yiqlisMrJdqfE5T58+TceOHav8fXv27CE+Pp6YmBhmzpyJsyJUqqoyc+ZMYmJi6NSpE3v37q1xvBURCVgQBEFwOYtVZqP5Evct2E7sC2u5b8F2NpkzsVhlj8Tz2GOPsWjRIo4fP87x48dZt25dmWPWrl17bXzhwoU89thjbo1JJGBBEATBpWS7wrbULCYu20NKei52RSUlPZcJy3bzS2q2S2bCACdPnqRLly7s2rWrwuPS09PJzc0lKSkJSZIYN24cX3/9dZnjVq1axbhx45AkiaSkJK5cuUJ6erpLYnVGbEMSBEEQquz/fXuYlAu5Tsf+M6oLszccdzo2+4djdGgWyBMr95UZax8ZyP8N71Cp65vNZkaNGsWHH36IyWQiISHB6XGbNm3i/PnzREVFXXstKiqK8+fPlzn2/PnzNG/evMxxTZs2rVRMVSUSsCAIguBSYQFGzBfznI6ZM/IICzDW6PyZmZncfffdfPXVV7Rv3x6AffvKJnRvJxKwIAiCUGUVzVQtVpm48ABS0svOkOMiAii2KXw2uWe1rx0UFER0dDRbt26lffv2mM1mRo4c6fTYTZs20axZM9LS0q69lpaWRrNmzcoc26xZM86dO3fT41xFJGBBEATBpQxaDbMGtWHCsrJtaZ+8vQ1GXc2WHxkMBv73v/8xePBg/P39GTNmTIUz4EaNGhEYGMiOHTvo0aMHy5YtY8aMGWWOGzFiBHPnzmXUqFHs3LmToKAgt91+BpGABUEQBBfTaTX0iglh8bjuzP7hGOaMPOIiAnjy9jb0iglBp635+l8/Pz9Wr17NoEGD8Pf3Z8SIERUeP3/+fMaPH09hYSFDhw5l6NChALz33nsATJkyhWHDhrFmzRpiYmLw9fVl6dKlNY6zIpKzvVA11b17d3X37rKffARBEIS668iRI7Rr167Sx8t2hWJZwUevpdBmx6jTuCT5epqzn4MkSXtUVa1S0yIxAxYEQRDcQqf9PeH6GUW6uVHd/ygiCIIgCHWQSMCCIAiC4AEiAQuCIAiCB4gELAiCIAgeIBKwIAiCIHiASMCCIAiCe9hlsBaAqji+2mveCcmd7QiXL19Op06diI+Pp1evXuzfv7/G8VZEJGBBEATB9awFcPx7+GAwvBzq+Jq6wfG6B1SmHWGrVq3YvHkzBw8e5MUXX2TSpElujUkkYEEQBMG17DKc3AwrR0PGQVDsjq8rRsGpLS6ZCYPr2xH26tWLxo0bA5CUlFSqfrQ7iJ3RgiAIQtWtfc6RVJ350yLY9KrzsY2vQtPO8N+JZcci4mHoa5W6vDvaEV5vyZIl18pVuotIwEKdI9sVrHYFk15Lkc2OQVs/ytsJQr3hHw6XUpyPXUoB/yY1Or272xFu3LiRJUuWsHXrVped0xmRgIVaV5MEarHKbEvNYvaG45gv5hEXHsCsQY4C774G8ddZEGpNRTNVawE0ae98htykPcjF8Mh31b60u9oRAhw4cIAJEyawdu1aQkJCqh1jZYh3LKFW1SSBynaFbalZTFy259prKem5TFi2m8XjutM/LkzMhAXBG2iNMOAFxzPfGw143jFeA+5qR3j27FnuvfdePv74Y9q0aVOjGCtDvFsJteb6BJqSnotdUa8l0F9Ss5HtSqnjbXaFi7lFHL5wlS3HMskrlpm94bjTc8/+4RjFsuJ07MYYLFYZRVWxWOUy1xQEwQW0OmjVF0avhIhOoNE5vo5e6XhdW/O5X0k7wtmzZ/PNN9/c9Pj58+czYcIEYmJiaN26dal2hCUtCV966SWys7OZOnUqCQkJdO9epeZGVVapdoSSJJ0G8gA7IN+s5ZJoRyg4Y7HK3LdgOynpuWXGOkQG8sH4W5m54leyC6xk5RdzxWIrdcyJV4fR5m9rsStl/87qNBJHXxnC9OW/EtPE/9qvW8L8rs2sxe1rQaiZqrYjxC6DvRj0PmArdMx8XZB8Pc0T7QgHqKqaVZWTC8L1THot5ot5TsfMGXmE+BtQVYht4k/SLcGE+hsJ8TcS5m8gxN+IxSoTFx7gNIHHRQSQnWfl2MU8Nhy5WCpJN2vkw5O3x+Jn1PHY8r3XXhe3rwXBzbS63xOuwc+zsXihuv9RRPBq6VcLWXMwg+8OXODd0V0qTKDFNoXPp/Qs91yyXWHWoDZMWFb27sqTt7chxN/AT0/3p1i2cybbQuql/Gu//hAbyp8/dH5XZvYPx+jZOkQkYEEQalVlE7AKfC9Jkgq8r6rqQjfGJNQBFa1kzrhaxJqD6Xx3MJ09Z3IAaNc0kHOXLTw5KLbUIqoST97eBqOu4gSo02roFRPC4nHdmf3DMcwZecRFBPDk7Y7byCXXN+q0tAkPoE14wLXvVVS1wtm3Sa+lWLZj1Gmr9fMQBEGoqsom4N6qqp6XJKkJsEGSpKOqqm65/gBJkiYBkwCio6NdHKbgTZw9S31yUCztmwbx7Jf72XoiG4C2EQE8fUcbhsU35ZYw/2vfe7MEWhFfg47+cWH0bB2Cj15Loc2OUXfzbUxFNnuFs+8LVwq5e9427u8exZjEaFqEiNtlgiC4V6UWYZX6Bkn6B5Cvquqb5R0jFmHVX7JdYaP5ktNZ7NwxXcjOLya3UGZYp6a0/i3pOjtHsaxUKYG6Iu5N5kynt68Xj+vOLWF+/Hud+drz4z6xoTzYowW3t2tyLTZRAERo6Kq8CKueqrVFWJIk+QEaVVXzfvvvO4CXqnIRof6w2pVytwIt2HSCzyf3xM9Y8V8r3XWJ62bHusrNbl/7GnS8N7YbGVeL+GzXOVbuOsuUT/YQHmjkoaQWjOvZkuRT2WIFtSAILlOZj+/hwFZJkvYDycB3qqqWbSMh1HuKomLSVbyS2Ufvvc9QS25ffz65J8deGcrnk3vSPy6sVAKNCDLx+O2x/PyXASwa1512TQOxygo/H8us9P5lQRAcZEXGYrOgqAoWmwVZ8e52hJs2bSIoKIiEhAQSEhJ46SX3zjVv+tFdVdWTQGe3RiF4vX3nrvCPbw4zd0zFK5kLbfZam9VWR2Vn3zqthkHtwxnUPpy8Ihsj39/h9DixgloQnLPYLOxM38m8ffNIvZJKTKMYpneZTmJEIr5631qPp6QdYY8ePRg2bBjr1q1z2myhT58+rF69ulZiEu8aQoUu5Rbx1Of7uWfeNs5fKSSnwMqsQc5LtFVmJXNd5GfU1dlZvyB4gqzI7EzfycyNMzHnmLGrdsw5Zmb8NIPkjGSXzITB9e0Ia5v3TlUEjyqW7Szddpo5Px7HZleZ0q810wfG4G/U0bqGK5nrmputoM4tstHI1+CByATBc15Pfp2jl486HXutz2vM2zfP6di8ffNoF9yO535+rsxY2+C2PJv4bKWu7652hNu3b6dz585ERkby5ptv0qFDh0rFUx0iAQulqKrKT0cv8fLqFE5nW7i9XRP+dmd7Wob+vi2nuluB6iqDVlNuAZCpA2JYsvUU5y5b+PvwDgT7iUQsCKG+oaReSXU6lpqTSqhPaI3O7652hF27duXMmTP4+/uzZs0a7rnnHo4fd77o1BVEAm6gnG2pOXPZwkvfprD5WCatw/z46NFE+rUJc/r9nljJ7CkVraDu2TqEk5n5rD6Qzs/Hs3jp7o4Mi49AkiRPhy0IblXRTNVisxDTKAZzjrnMWEzjGIrtxSwdsrTa13ZXO8LAwMBr/z1s2DCmTp1KVlYWoaE1+8BQnvr9zik45ayQxszbYogIMnEkPZe/3dmOh3u1RF9PZ7TVUdGsf8bAWG5vF85fvjzAtE/3MrhDOC/f05EmASZPhy0IHmHQGpjeZTozfirb8m9awjQM2prdKXJXO8KMjAzCw8ORJInk5GQURXFrT2DxDtvAlNcScMone0m/WsQPs/oxoc8tIvk6odNq8DPq0Ggk/Iy6Urfc2zUN5H9Te/HskLZsNGcy6O0tfLknzelWB0Go73QaHYkRicwZOIe2wW3RSTraBrdlzsA5JEYkotN4ZzvCL7/8ko4dO9K5c2dmzpzJypUr3Xo3q8qVsCpDVMLyXjdrCViZQhpCxU5k5vPslwfYfSaHfm3CePXeeJo18hGVtIQ6r6qVsGRFxmq3YtKZKJKLMGgNLkm+nuaqSljiX38Dc7OWgGJLTc21DvPn88k9+cfw9uw6fZl75m7lwpVCNpovcd+C7cS+sJb7FmxnkzkTi9U12zEEwRvpNDp89b5oJA2+et96kXxdSSTgBiavyNFT15mSQhpCzWk0EuP/0Ir1T/TlmSFx7D2TIyppCYJQikjADciK5LN8sPUk0wbGOB2vr4U0PKl5sC93xUcyf9MJp+OzfzhGsSwSsCA0ROLdtgGQ7Qr/+OYwf/3qICnpefSNDWXxuO50iAxEp5HoEBnI4nHd620hDU8zGcRtf0EQyhI35Ou5qxYb01fs5efjWfy5dyv+OrQtOq2mQRXS8LSbVdLy9vrZgiC4h3jHrcdOZubzx/nb2HEym9f/FM+Ld7W/lmQr2lIjuFZJJS1npg+IQacRRTsEoSES77r11NbjWdwzbxtXCm0sn5DEyFujPR1Sg3V9Ja3rb/u/91BXmgaZmLp8L/nFYjW0UP8osoxisaAqCorFgiJ7dzvCq1evMnz4cDp37kyHDh1YurT61boqQyTgekZVVZZtP83DS5NpGuTDqml/ILFVsOsvZJfBWgCq4vhqFwmkIs56Ed/eLpyTWQVsOpbJqIXbycwr9nSYguAyisVCweYtnB7zIEfjO3F6zIMUbNmCYrF4JJ6SdoTHjx/n+PHjrFtXtq39vHnzaN++Pfv372fTpk089dRTWK1Wt8UkEnAdJtsVLFYZRVWxWGWKZTsvfZvC31cdZkBcGP+d2ovmwZXsu1mVhGotgOPfwweD4eVQx9fUDY7XhXI5u+1/b9coFj/cnROXCvjTgl84nSV+hkLdp8gyBdt3kDZtGsVHj4LdTvHRo6RNnUbBjh0umQmD69sRSpJEXl4eqqqSn59PcHAwOp371meIlR91lLN6ztMHxnB3QiQBPjoev60N2so+W7QWwMnNsOlVuJQCTdrDgBegVV8w+JU+1i47jl05+vfXMg7CilEweiXEDAKtk79WdhnsxaD3AVshaI3Oj2uABsQ14dOJPXj0w13c994vLB2fSHxUkKfDEoQKZbz6KsVHnLcjjHzzDTLnzHE6ljlnLqb27Tn3zF/KjBnbtSXi+ecrdX13tCOcPn06I0aMIDIykry8PD777DM0GvfNU8U7YB10fT3nEinpuUxdvpf3x3Zj5sDYyiffihLqyE9A7wspX0NBFhRkwv0fOhK1MxtfheiejiSrM/7+elUSfAPVJboxXz7Wi3FLkhm1cDvvje1Gn1jnnagEwdvpQkMpLqeNX/Hx4+jCavZ3213tCNevX09CQgI//fQTJ06cYNCgQfTp06dUlyRXEgm4DrLaFWZvcP6X+90fj9M7JrTyq5rtxeUn1M3/diThY+vBLwz8QiGgqSOJOnMpBYz+8M9IiOwCzROh82jIOQUrx/x+XGVmzA1Q6zB/vprai4c/SObRD3fx5v2duTuhbMs0QfAGFc1UFYsFY2ys4/bzDYyxsajFxbT4eFm1r+2udoRLly7lueeeQ5IkYmJiaNWqFUePHiUxMbHasVZEPAOug1xSz/nqedj8BuhMFSfURs3h6WPw2DYYtwrkIscM1pkm7aEoFxInOp4l73wPjH6w6V/Oj9/4quMDgHBNeKCJz6f0pFuLxjy+ch+Lfz7p6ZAEoeoMBsJmlm33BxA2YzoYXNOOcNmyZXz66afExcWxb98+p78aNWpE06ZNr7UjVFWVZcuWcffdd5c5b3R0ND/++CMAFy9exGw2c8stt9Qo1oqIBFwHlRR2cKZz8yBk2eZ8QZXdBkdWw/IH4J2OsPEVyL9YcUK1FZZ+TWt03D52ZsDzYAyEwf+ECRvguXMQ1LziBK/3qcSfuGEJNOn58JFEhsVH8Mp3R/j3uiPYZHupBXeifrTgzTQ6HX5JSUTNn4exXTvQ6TC2a0fU/Hn4JSWhccHCJne0I3zxxRf55ZdfiI+P57bbbuP1118nNDS0xrGWR7QjrIMuXi1i95nLTPv011Kv+xq0JD/TE7/z25Cuf97a/6+Or8vugSunHbeREx6ELg85EmTqBsct4RuVd4vYWgCntjhmsNee6T7v/JmutcCxSjrjYNnzR3SCh78Fn0Y1+4HUU3ZF5V9rjnBXp6ZcuFrE3J9Sry24mzWoDb1iQvA1iNv3Qu2pajtCRZbBakUymVCLisBgcEny9TRXtSMUCbiOOZttYeySncx9sCsZVwt554fjmDPyiIsI4N2RnbglZxvS9QuqStz/IchWMAVBzO2lk2pVEmqJyq5qtsvlJ/j7P4SLhyHzKNz2Dwh13iSiIbPZ7fxw5BKPfbK3zNjicd3pHxcmqpgJtaaqCbi+clUCrvsfRRqQS7lFPLRkJ7lFNkx6DQPimtCrdei1es6+UjHS/8pZUPXz2/DoOucJ1eDnmOm26lv5bUJa3e/jFa1k1uoc5x29smyCb9kbLp+CExvh6Bro/gj0exb8m4htS7+x2VXm/JjqdGz2D8fo2Vo00BCEuqrhvaPVUVctNsYuSSYrv5jlE3oQ28TxDLjkzdfPqANVU/3nrZVNqNVRUYLvM8txK3zz67B7KaR8AxN+cMyMxbYl1yy4EwTBK4mPznWAxSrzyIfJnMoqYOHY7nSJblz2oIuHITe9aguqapNW50ieksbx9frZrH8TuPMtmLbTsYDr/B7HvuSMg6DYf9+2dGpLgyt5WdGCu5JOSoIg1E0iAXu5YtnO5I/3sO/cFd4dnUDv2BtW5Cl22PYfWNgfDnwO/Z5zfqIBzztmnd4sNBba3glb33Y+3gC3LVXUSWnGwFiMOvFPWBDqKvGv14vZFZVZn+3n5+NZvHZvJ4Z0bFr6gMsn4cM7YcPfIfYO6DoWWg9wPG+N6AQanePr6JWO27d14Rmq3kdsW7pOeZ2UFjzYlYhAIxvNmZ4OURCEaqr0O7IkSVpgN3BeVdW73BeSAI6uRn/7+iDfHUznhWHteODW5tcPwp4PYf0LoNHCH9+HTiNB+q38ZFUXVHkTW6HjdrmzbUslt9Eb2HPgkk5KPVuHXFtwp5Ukxn+YzK9nr7ByUpLzxxKC4GGKXcEuq+gMGmSrglYnoanhosHTp09z1113cejQoSp93wsvvMCyZcvIyckhPz+/3OP+9a9/sWTJErRaLe+++y6DBw+uUbwVqcpP4nHgiLsCEUp7fZ2ZVfsusPzRbkxMCv+9qIblCnw2FlY/AVHdYep26Dzq9+QLFT9v9XYVFfroMwtyTtdqON7ixk5KJoOWuWO6Eh5oYsJHuzmTLbooCd7FVmznzKFsvnpjDwumbeKrN/Zw5nA2tmLPrFsYPnw4ycnJFR6TkpLCypUrOXz4MOvWrWPq1KnY7e6Lt1IJWJKkKOBOYLHbIhGueW/zCZZtP83GmbfSS9lbuu3fqY3Q50m4czaM/RqCom56vjrl+m1L199GH7UCQtvAkjtg+3xPR+kVQv2NLH3kVuyqyiNLd5FT4L6+pYJQFYpdIe3oZdYsOEhWWj6qopKVls+a+QdJM19GcVElt8q2IwRISkqiadOmFR6zatUqRo0ahdFopFWrVsTExNw0addEZadG7wB/AZwvxwQkSZoETAJHPU2helYkn+W1tUeZP7oTTbJ2li6qkXEQvhgPoz6FruPAjW2yPKq8bUuKzfGMe/1f4cpZx4ppTcPehtM6zJ+FY7vz0OKdTPp4Nx//uQcmsTVJqAU/f36MrHPOb+UO+nMHkr895XQs+dtThDUPYMMHZdd6hDb3p88Dzhcd3qgq7QgbNapctb3z58+TlJR07ffltS10lZsmYEmS7gIuqaq6R5Kk/uUdp6rqQmAhOCphuSzCek62K1jtCia9ltxCGxeuFDK4QzhD2wUjfVBOUY1Nr8Et/evWreWqcrYvWauD+z+C71+EHfPg6jm4dxEYfD0XpxdIbBXMWw90ZsaKX3nmywP8Z2QCmsq2oxQEN/ANNJB9wfljkcvnC/ANrNmODHe1I6xtlXkH/wMwQpKkYYAJCJQk6RNVVR9yb2j1n8Uqsy01i9kbjl+r8Tt9YAyT+7VG0tegqMZvZEXGardi0pkokoswaA3oNHU8aWu0MORVaBQN656D5ffD6BWO1+viojMXGd45krScQl5fd5Soxj48O6Stp0MS6rmKZqq2YjshkX5kpZWdIQc388MuK/zxqa7VvnZV2xFWdgbcrFkzzp07d+335bUtdJWbvkupqvpX4K8Av82AnxbJt+Zku8K21CwmLttz7bWU9FymLt/L4nHdue0WH6QarAa22CzsTN/JvH3zSL2SSkyjGKZ3mU5iRCK++nowY0ya4kjCAeGQ+iNsfavBV82a0u8WzuVYWLDpBM0b+zKmh3gUJHiGVieROKIVa+aXff9KHN4Kra5md2hK2hEOHjwYf39/xowZ45IZ8IgRIxgzZgyzZs3iwoULHD9+3G29gEHsA/YYq11h9objTsdWrN+MuudD6PNU2UGDH/KfFmGRJBRVwWKzICulq0PJiszO9J3M3DgTc44Zu2rHnGNmxk8zSM5Idnq8xWYp93xeK/YOyEuHL8eLqlmAJEm8NKIDA+LCeHHVITaaL3k6JKGB0mg1RMUFM2xqPKHN/dFoJEKb+zNsajxRccE13ooEVW9H+Je//IWoqCgsFgtRUVH84x//AOCbb77h73//OwAdOnTggQceoH379gwZMoR58+ah1bpvTYXohuQhiqoS+8Ja7Erpn38rKZ0Vhn8SHuSDNOVnOJf8exODyG5YRn/CzuzDTme2eo2e07mnCfMNY8L6CZhzzGWu2za4LR8M/gCT1oReq6/bM+WbtTosr/lEPVdQLPPA+9s5lVXA55N70rFZkKdDEuqJKrcjLNkHrNcg21yzD9gbuKobUt3/SdRRzmr8xkhpfGZ4GZPGTuH9K8E3xLEa+NF18GIm8sOr2Jl92OnMduv5rSw8sJB7v7mXAH0AqVecd9BJzUnFR+dDj0978NnRz9h2flulZ8peR1TNcsrPqOOD8bfSyEfPtOV7ySmwYrHKKKqKxSoju2gLiCDcjEarQW/UImkk9EZtvUi+riR+Gh5y8lI+Uwe0vvb7OOksKw2vAHDkjhUYIjs6Bq4rqmFFZd6+eU7Pt+jgIka1HcVrfV7DIluIaeS8t25M4xjybfk81O4hBkQP4P0D7zs9bt6+eVjtXr6vtKRqljOebj7hYeGBJj58NJE5o7vwy4ks7luwndgX1nLfgu1sMmdisXr5hytBaABEAvaA4xfzeHjpLmKb+LN4XHeGN8lkheEV0OhJHfYZnbv1KNPj1WKzYNKZKpzZNjY15s5b7sRH58P0LtOdHjctYRr+en9mdZ9FqE9ohecz6Uw1+4O6W0VVs/o+DVpD7cbjZW4J9SM9t4hpn/5KSnoudkUlJT2XCct280tqtpgJC9XijseWdYkr//wiAdcyi1Vm6vK9ADT2NdA/4Bzvqq8ROPg5fJ7bT1JiIkjWa7d/rXYry48sZ+hXQ0kvSK9wZlskFwGg0+hIjEhkzsA5tA1ui07S0Ta4LXMGziExIvHaVqQiuajC82UXZnPJ8vtCHq9brFVe1ayRn0BgM/hlrmfj8zCrXeE/Pzhf6Df7h2MUyyIBC1VjMpnIzs5usElYVVWys7MxmVwzOWlYmyU9zNFg4RCpmfl8/GgPmlzZD5+Pw/LQl+y0ZjJv/SOlFkK1C27HpO8ncTL3JLdG3IqqqkzvMp0ZP80oc+5pCdMwXDfj89X70rtZbxIjEsvdB2zQGso93+ROk/nC/AUfpXzEzK4zuSfmHpLTk71vsZbTqlkGWP0U/PoRhNwC7e/2XHweZNJrMV/MczpmzsjDR1TMEqooKiqKtLQ0MjMbbhcuk8lEVJRrSgCLVdC16PC5y0T4SwQHBSEV50LyQmSfRvwcEcPMjY+XOf7Nfm+SU5RDdGA0PZv2RJIkLDYLyRm/JcKcVGIaxzAtYVq1E2FF58suzOa1Xa/RLrgdsY1jeXrz02W+f87AOfRu1tv7CnzYiuCju+DiYXh0PTTt5OmIap3FKnPfgu2kpOeWGesQGcjnk3viZ/Sy/2+CUEdVZxW0SMC1pMiSByc3Y9r62u8FI/o8jaXNIMatfbjcLUMfDfmoTGJ1dYWrm53vavFV/rz+z1WK0SvkXZ2fFUYAACAASURBVIRFAxxdoSZuBP8wT0dUq2S7wiZzJhOWlf23uGhcdwbEhZVZayAIQvWIbUheqthqRTq5GdOXD5YuGPHFw5h0PlVeCKXT6PDV+6KRNPjqfWs8+7zZ+QIMFW9r8trFWgHhjsYVBVnw2UMgF3s6olql02roFRPC4nHd6RAZiE4j0SEykLljuhDbxF8kX0HwMPEv0M1UVcVWXIhx62tOx4vyMyq1sMqTbrZYyxtiLFdkAvxxAZzbAatnQQNbPOJr0NE/LozPJ/fk2CtD+XxyT3z0Woa9+zMbj4pKWYLgSSIBu9nynWfx9QtwWjCiWILNKZ8xsdNEp99748IqTylZrOXM1M5TvSLGCnX4I/R7FvZ9AsmLHCUqrQWgKo6v9bxkpU6rwc+oQ6OR8DPq+ENMKNHBvjzz5X4y8xrWXQFB8CYiAbvRofNXeenbFAryc8sUjMjRaJgQEc7/pa4gPqTjTbcMeVJ525re7PcmTf2bcrX4qqdDvLl+z0H8/dCsKxxb5yhh+XKo42vqBkcibiBMei3vju5CXpHM01/sR1Ea1l0BQfAWYhGWm+QW2bjr3a1YZYWf75XR2y3wxXgAzuh0TI0II0Or49W4sQzu8RSyROmFUOjQ2GQkkwm1qAgMBjS635OxIstgtZY77g43LtbKKszioTUPEeITwpLBSwg2Bbv1+jVmtcDx9df+P5QyeqVjO1MDamP48fbTvLjqMH+/qz2P9m7l6XAEoU4Ti7C8hKqqPPvlAc5fKeTTpHPov3rEUSBi9Er2RrbnwcgI8nQGlsRPY3D36aDVlVoIZbJB4ZatnB7zIEfjO3F6zIMUbNmCYrEAoFgsFGzeUu64u9y4WCs6MJo3+r3BubxzTPx+IjlFOd5XrKMUFX5+y/nQxlfB3rBuxz6U1ILb2zXhtbVHOeJkq5IgCO4lErAbfLLjDDHhvux+oRct+4/FMmMPclBz1utUJvgU07hxK5YP/4KELhPLdOtRZJmC7TtImzaN4qNHwW6n+OhR0qZOo2DHDuxWa4Xjily7Ca9H0x7MGTiHS5ZLXLJcYkvaFsatHUfXj7sybu04tp7fisXm3g8GlSaaN5QiSRKv/6kTQb56Zq74lSKb3dMhCUKDIhKwix1Iu0KnaF/iYy4w8ccJdP2kG+M2zmDT+S00C4zmH4l/5393fEZU41iUYmvZhGm1kjlnjtNzZ86ZW6nx2tYzsifLhi7j1NVTPL7xce/trCSaN5QR4m/k7Qc6c/xSPv/87oinwxGEBkUk4BqS7cq1Vm8FxTJFNhsZxQd5assTpRLR81ufJ9rQhAFpgZx98KHSt44LCsj/ZTsX/vo8ktFI8XHn9XuLjx9H4+tb4bjkohqlVRXuG87ig4udjnlNZ6WKmjcMeN4x3gD1iQ1jQu9WfLzjDBtSLno6HEFoMEQCrgGLVWaj+dK1Vm/3v7edNk19WHRoQZljx7cdiy15r9Nbx/lbfqZw1y7yNmzAfvkyxthYp9czxsaiFhZWOK7k5wOOW9mKxYKqKCgWi9tvTd+sU5NXFOsor3nDAx87Xm9AC7Bu9MyQONo3DeQvX+7nYq4X7+sWhHpEJOBqku0K21KzmLhsT6lWbwEGX6eJaGSre7k6z3nv3ayFC2n88Dhit21F07gxYTPLNkcACJsxHVWvL3c8dPIkLn+ynOITJyjYtLlWF2nVmWIdJc0bHl0HL2bCg184qpJln/B0ZB5l1Dm2JhXa7Dz1udiaJAi1QSTgarLaFWZvKH0rOFZKoyjvvNNE1LhRRIW3jrWBgWiMRjQ6HX5JSUTNn4exXTvQ6TC2a0fU/Hn4JSWhNRjKH+/dB9/u3Sg+doy06dNrdZFWRcU6vKWgyDVanSMRSxrQGWHvR7BqWr0vyHEzMU38+ftdHdiamsWSrac8HY4g1HsN955bDd3Y6i2QAhbq30Lek8/E+Mk8vWVWqeOv5GRgjI11JMUbGGNjUYuKkHwdDQ00vr749e2LX1KS032+FY37dOzI6TEPOo05c85c/JKSwA37ha8v1nF9Z6WJ8RO9pqCIUz6NYdgb8Pk42D4Xej/h6Yg8anRiczYfu8S/1x+lZ+sQOjYL8nRIglBviRlwNRXZ7MSFBwAgofCOfh5RUhYP78qkqX84b/d7+1rVqPvPRCCv/JrQyZOdnitsxnQwlJ4hanQ6NL6+SBoNGl/fMkU2yhuXTCaPLdIq6UH80ZCP2DN2DwsHLeRYzjG+Tv3abdd0iXYjoO1dsOlfDf5WtCRJvHZvJ4L9DMxc+SsWa8O+KyAI7iQScDUZtBqeGORYDPWM72r69x9E4YwjfPXMP2mua0Lfpr1Z2u89vjPfyf2fpmHZvRu/P/Qq99ayq6pYqUVFFS/iKnLvs9jri3U0MjbCnGPmrd1vkZrjfIGWV5AkGPYmaA3w7eMNrmHDjRr7GZj9QAKnsgp4eXU5+6YFQagxUYqyOuwyqr0Y9D7IBTloC/PJ3XuEywvep/j4cYyxsYTNmI4+MpLTD40l+KGHCJsxHUmnK1VCUrHaUCQtOoMW2aqg1UloatgiTpFlCrZsIW3qtDJjzWbPxufW7uhDQ2t0jarILszm3m/uJdw3nOXDlqPX6mvt2lW2eymsfgJGzIGu4zwdjce9tvYoy7af5qvHehEd4otJr6XIZseg1YhWhoJwg+qUohQJuKqsBXByM2x6FS6loNzxOgX5LUmbPrPUYZKvLy1WrETbvAV6H32ZBGsrtpN29DLJ354i+0IBIZF+JI5oRVRcMHqjtkYhKhYLBTt2kDln7u8fCKZPQxcRwflZTxG9eBG6pk1rrZb0T2d/4vGNjzMxfiIzu868+Td4iqLAR3fBxUMwLRkCIjwdkUdZZYUTmfmcyMxn/sYTmC/mERcewKxBbegVE4KvwUuf6wuCB4gE7G52GY5/DytHX3tJmX6I0xNmllpcJfn6ErVsOZcKA9i1/nyZBKvVSZw5lM2aBQfLXGLY1HhadAhxyUz4xgRrO32as5MmEzV3DnJ6eukEPXOG41b4bwvBXO3FbS/yzYlv+GjIRyQ0SXDLNVwiKxUW9IK4IfDAMk9H41GyXeGno5eY9PGeMmOLx3Wnf1yYmAkLwm9EMwZ3sxc7Zr7XkYIjyyx6Ch7/CBct/qxdYiYrLR9VUclKy2fN/IOkmS9jlxWSv3W+zSP521PYZRXFrmArtqOqKrZiO4pdqVKozhZpGWNiaPnpcmxnzpA2rXa3KT1767O0DGjJwcyDFNgKvLRZAxAaA/2fhZRVcGS1p6PxKKtd4Z0fnC/om/3DMYrlqv2dFAShNJGAq+KGYv7r/Hy5fCWjzKKngPtHsvv7C05PkfztKXQGLdkXnPefzc1yLJI6cyibr97Yw4Jpm/jqjT2cOZyNrbjmxfK1gYFkvb/Q6Zg7a0n7G/z5YMgHhPmG8fDah72zWUOJXjOhWSJkp0JxPqiK49FDA9snfONWu+uZM/Lw0dfsUYkgNHQ3TcCSJJkkSUqWJGm/JEmHJUn6f7URmDdSrL8X8z9gNPC30GA2H/uC0CmTSh1nahJSboK9fL4A2aoQEunndDzp7lacO3KZNQsOOp09V3UmfCNPbVOSFZkDmQd4Zssz3tusoYRWD2NWQKNoWDoEXg6FDwZD6gZHIm4grt9qd6O4iAAKRfckQaiRysyAi4GBqqp2BhKAIZIkJbk3LO8jKzJ5qooyeQv5Tx5i1x0v0lwykrhgHvqmkTT7zzvXthcVXb5aboINbuYHEiSOcN4APfbWCHatrvj2dE14apuS1W5l3r55Tse8pllDCbsM55Lhy0ccZSoVu+PrilFwakuDmQkbtBpmDWrjdOzJ29tg1IkbaIJQEzf9F6Q65P/2W/1vvxrURkmLzcLmc1v48/fj6fpxV8ZvfpIWIe35MH4eIa/9gik+Hp9+txH96Qpsc77jUPIVug1t4fRcicNbodVKRMUFM2xqPKHN/dFoJEKb+zNsajxGX12Fs2edvoZvegZDhbWmbywI4ip1ollDCSfP+q/Z+KpjvAHQaTX0iglh8bjudIgMRKeR6BAZyLwxXUhqHSIWYAlCDVVqH4EkSVpgDxADzFNVdaeTYyYBkwCio6NdGaNHyYrMzvSdPLHp8Wuvnc07S2tNFFl2f3YvOUX2hUOERPrRbWhLwlsGse/Hc8T3j2LY1HiSvz3F5fMFBDfzI3G4YxW0RqtBo4UWHUKIigtGp9cg2xzblEpuT2el5ZeJJbiZH7JNqdE2petrTV+/Cjr0sSkuLQhyo5JmDeYcc5mxkmYNvnr3rMCushue9ZdyKcUx3kD4GnT0jwujZ+sQfPRa8optLP75FCcy85l5m/PZsSAIlVOpd1tVVe1AgiRJjYD/SZLUUVXVQzccsxBYCI5tSC6P1EOc3Tp9MuphlAx/1n1w7NprWWn5rF90iKFT4rl9fHu0Oo3TBHv99qKSRAxcl1QVEke0Ys38sluUug1pQf6VIoJCfbDLKjqDploFPG6sJW2/fJmc5Z+iaxKOb0Lnyv9wqqCkWcOMn8rOvr2uWYPtt2f9GWX/H9CkvWPc4PwRQ32ku67wRpCPgbScQr47kM6dnSJpHebv4egEoe6q0j0kVVWvABuBIe4Jx/vceOs0Mlvl7rgHyl3lvOu7Uyh2x+cPjVaD3qhF0kjojdpKJUmNVuP89vRj8TSK8MVaZOf0gawar5AutU3Jx4erq1aR/rcXUNy0Cvr6Zg0lNbLbBrflzX5v0i28m3c1a9AaYcALzscGPO8Yb8D+OqwtRr2G/1t1GHfUERCEhuKmhTgkSQoDbKqqXpEkyQf4HnhdVdVyN0nWl0IcdkXFUmjBbssjMDCMKzkZ5H7xX1pMnM6CGZtRnfRM1Wgkpsztj6SRanRtxa44ZrnX3562KZw9nM36RYfLHF/TAh75mzdzbvIUQmdMJ2xa2TKWriIrMla7FZPORIGtgGWHl1FsL2ZW91k3/+baZC1wLLja6Kh4RpP20GcWtL4NTIGejs7jlm0/zd9XHWbO6C4M7xzp6XAEwePcVYijKbBRkqQDwC5gQ0XJtz65cjkX5ZdfuDz+Mczxnbny5+k0bhFLUZGtwlXOsq3mBQqczZ4lSWLP2jNOj6/pCmn/fv0IuuceUFTseXmoioJisbi8MMf1zRoCDAFcKrzEJ0c+IT0/3aXXqTGDH8QMgkfXwYuZMH41XDwMO5yv5G5oHuzRgo7NAnl5dQp5RTZPhyMIdVJlVkEfUFW1i6qqnVRV7aiq6ku1EZinXc0vhD07uTBjZqmKUcdfeJ3Dm9IqXuWsq9nstzw6g8atK6TDX3geY0xrzowdx9H4Tpwe8yAFW7agWNxXKGNKpykAvHfgPbddo9q0OkciljRgCoLLp2DbHMi/5OnIPE6rkXjlnngy84vLrZYlCELFxD4CJ1RVBauVqwsWlHrd4tOEXzs/zsG1ZqLiGjndRlSyytkdKirgUdOZtyLLWJJ3cf7JWbVaorKpf1NGxo1kVeoqTl11vv/Zawx4HuQi+PltT0fiFRKaN2J0YjQf/nKalAu5ng5HEOockYCd+Gb/BQKD/EtVjLKYQvk1YSaqpCF+59sYffW06BDCvU93Y8rc/tz7dDdadAipcSejimh1UrkFPGo887ZayZwzx+mQO0tUAkyIn4BBa2D+vvluu4ZLhLSGLg/C7iVw5Zyno/EKfxkcR5CPnhdXHUJxsiZCEITyiQR8g5wCKyu+WYM158K1ilGFphB+TXgcRaOny/53CY4KRC0qqtYq55oob4X04IkdaBbbuEbX91SJSoAQnxDGth/LutPrOJJ9xG3XcYl+zzq+bn7ds3F4iUa+Bv46tC17zuTw5Z40T4cjCHWKSMA42q5ZrDKKqiLLVuYPULGZGtPqq69ouWkryovvIfn6krB/Dv4FF9xaMepm9EZtqZn3PU92Ift8Adu/dl5lqrI8VaKyxMMdHibQEMicX53Pwr1GUBTcOgH2fQpZ4tknwJ+6RnFry8b8a+0Rcgq8qKSoIHi5Bp+ALVaZjeZL3LdgO53/3/eo+VfI9h3Gd++ksGD6Jr5ZfIJGLcK4b0YcYbFhRM2f59aKUZVx/czb6KtHVVVMfgaKLbZqty/0VInKEoGGQCZ1mkT7kPbkFud6b7tCgN6zQGdybFES0GgkXr6nI7lFMv9ef/Tm3yAIAlCJfcDVUVf2Act2hY3mS0xc5mg4/n8Dw+kbEs66RWVnNsOmxhPdJhBJp/Vo8nXGWiRz9nA2e9aeIftCASGRfiSOcJS9rMozacVioWDHjtIlKidPwq9PH7R+7q/8VGArYGvaVhYdXETqlVRiGsUwvct0EiMSvadMZYmfXoEtb8Dkn6FpJ09H4xVeWZ3Ckm2n+O9jvega3djT4QhCrXLXPuB6y2pXmL3BkWwlFEZ2b8nutc73oyZ/ewpF0nld8lXsCufNOaxfdLjG7QtLSlS2XP4JbQ/sp8WHSyk2HyPv++/d+CdwkBWZ5PRknt7ytPe3KwToOR1MjRyJWADgiUFtaBJg5MWvDyHXsG2mIDQEDToBm7Qqn4zryMlXh3Lo2URMjYLd24nIDeyySvK3rmtfWKpEZWAgeZs3kb1wEari3jfUOtWuEMCnEfR+Ao6vh7M7PB2NV/A36vj7XR04fCGXT3Y4LxgjCMLvvC+j1BZrAdLxDQR/NhzNK6Hs/2IEV65ku73Clau5sziHJEmETpiA9dQp8n78sdrnqYw61a6wROIk8A+HH18GURMZgGHxEfSJDeWt749xKde9C/cEoa5rmAnYLsPJzUgrR0PGQfJVhbcKc8n7/Cu6Dan9Clc14c7iHAABd9yBPiqK7MWL3Vp4v6RdoTMl7Qq9jsEP+j4DZ7bCyY2ejsYrSJLES3d3pFhW+OcaL99SJgge1kATcOmG628GN+ZP62DXNjsBIUaGTqndClc14dbiHICk0xH86CMU7T+AZdeuGp2rIiXtCp3xunaF1+v6MITGweWTjgYOquL4aveyZ9a1qFWoH1P6t2bVvgv8kprl6XAEwWs1yFXQst2GNT8DU0AkefkX+Om9f+K/TuFYm1H0DV1O+38sRLFTbh9fb2MrtpNmvkzyt6e4fL6A4GZ+dBvSgugOIRhMNV80phQVkTrwNkwdOxC9cKELInbOYrOQnJHMvH3zSM1JJaZxDJM7TaZXZC/vWwV9vdwLjufAW9/+vXPSgBegVd8G1Tf4ekU2O3fP3cbIW6MYeWs0PgYtRTY7hut6CwtCfVKdVdDetaS3FlhsFnam73S8yV9JpYetOX/+qpg9nZ8m2rSPjs2OISlFaI2ON053lpZ0lZLiHFFxwej0GooLZfb/eA5bsUK7Xk1rfH6NyUTwI+NRLYXY8/LQ+Pk5CnMYDC5dFe6r96V3s94kRiRi0pnItGSy+uRq+kX1c9k1XM4uw4V98OUjv7+WcRBWjILRKx0dlbQN7p8ZJr2WTyYksvPkZe5/bzvmi3nEhQcwa1AbesWE4GtoeD8TQbhRg/ooKisyO9N3MnPjTMw5ZlS7zLDlpzkUOxaNj0S/O2xIA+tmw/XSxTl0nD6Yxa8bzrrsuW3j0aMxtol1e6ek69sVmnPMvLP3HX44+4NLr+FSNzzOKGXjq47xBki2K+w7d4XpK34lJT0Xu6KSkp7LhGW7+SU1W2xTEgQaWAK22q1cys9g0/B17B+7j63D1uPT40UK/ZqxN34duoGTHbcN6/iMRZIkEm5rTk56AedSLtf4fIosY9mZXOudkno3602UfxQrj650y/ldQu/juO3szKUUx3gDdP0e+xvN/uEYxbJIwILQoBKw0Qa3pTXi6qPTMcd35vCMNzieE0rzHiYO+mzHZAysN8/sYrqH4xtkYN+PLuja46FOSRpJw8i4key9tBfzZbNbrlFjtkLHM19nmrR3jDdAJr0W88U8p2PmjDx89N7/aEcQ3K3BJGBFlrFs307mTMcszqrx5YB/f/zz0+jR9CpPJMzwzq0u1aTVaUi4vTnhLQMpLpSrXyMaz3ZK+mPsHzFqjaw0e+ksWGt0LLhyZkDdfJzhCkU2O3HhAU7H4iICKLTZazkiQfA+DSYBO2ZxcwFQkUhpOxa71kT7lA8peH8RI5oP9d6tLtXUoXczgiP9+PqtvSyYtomv3tjDmcPZ2Iqr9ubnyU5JQcYghrYayncnvyPX6oVN37U6x2OL0SshohNodI6voz6tF48zqsug1TBrUBunY0/e3gajruG89QhCeRrMvwLVYKTpewtpe/gwLTb/wi1j7yTu/Hf4W9IdjQf8AtFp6s+bpWJXOH8sh+8X17xGtKc7JY1uO5pCuZBvUr9x63WqzeDnWO386Dr4WyaM/MSxH7iePM6oDp1WQ6+YEBaP606HyEB0GokOkYHMe7ArPVuHiK1IgkADScC2YjtnD1/mm8WnWDBjM6sXHSMk0pfEf09D8vWtlX63tc2VNaI1Oh1+SUlEzZ+HsV070OkwtmtXa60Z24e0p1NYJz4zf4aieuniHa3OkXA1Gsd+4P9OgIJsT0flUb4GHf3jwvh8ck+OvTKUT/7cg6Ppufxw5KKnQxMEr1DvE7BiV0g7epk1Cw6WmgmuX5zCJYs/weMfqZVZXG1zdY3oUp2S9u8jas676Jo2ReNbOwUyRsWN4nTuaXak14HGBz0eA7kI9nzg6Ug8TqfV4GfUodFIBPno+enoJd5Yb8YqVkELQv1PwBXNBHdtuECjR/5cK7O42uaOGtElnZKQJM48NJasd/5T0zArbXDLwQSbgr17S1KJJm2h9W2QvBhkL+vi5EEajcRfhrQlLaeQFclnPR2OIHhcvU/AN5sJ6v1MtTaLq03urBEtaTQEDb+L/K1bkS/XfJ9xZRi0BkbGjaRN4zbkW/NRVAWLzeJ9fYJL9JwK+Rlw+H+ejsSr9I0NpUerYOb8lEpBsZf+vxOEWlLvE7Ct2F7nWgy6gkarISoumGFT3dNYImjECLDbyV2z1kUR39xD7R4iplEM49eNp+vHXRm3dhxbz2/FYnNtNS6XaH2bo0nDjnmiVeF1JMkxC87KL2bpNud3pgShoaj3CfjIiSN0GxLldMxbWwy6SkmN6Huf7sbkuf0ZOjmeiFZBLqlvbYyNxdiuHVe/qZ2VybIis+fiHp7Z8gzmHDN21Y45x8yMn2aQnJHsfTNhSYKkKZC+H85u93Q0XqVbi8bc3i6c9zefJKdA3KIXGq56nYAzc4vY+sVkAhsbGTyxQ51pMehKJTWi8y8X8fHftnNke7rLzh00fDhFBw5QfMr9Mxmr3cq8ffOcjs3bNw+r3QvfyDuNAp/GsGO+pyPxOs8MjiPfKvPe5hOeDkUQPKZeZ5///fcdAg915ut3D9E41Mi9T3djytz+3Pt0N1p0CKkTnY5cJTDUhyYtAjix55LrznnnnSBJ5H672mXnLI9JZyL1SqrTsdScVEw691XjqjaDL3R7BI5+BzmnPR2NV4mLCOCPCc348JfTZFytX1sABaGybpqAJUlqLknSRkmSUiRJOixJ0uO1EVhNnci4TNG+L7EE3Umg8SrB0UHXugXpjdp6P/N1pnXXJlw6k0dutmvqE+vDm+DXM4mr337rsq5L5SmSi4hpFON0LKZxjPeWEU2cCJIGdrqvj3Jd9eSgNiiqyn9+dF7mVBDqu8pkIRl4SlXV9kASME2SpHKqz3uP7z5/ntC0u5B1RgbN6Ick1d9nvZXVumsTAE7szXTZOQPv+SNBd92FkpeHqigoFotbuiMZtAamd5nudGxawjTvLSMaGAkd/gh7l0GRF5bS9KDmwb6MSYzm893nOJXlfKeCINRnN03Aqqqmq6q697f/zgOOAM3cHVhN7D5yEu2BY+SE9iYmTiKkeZCnQ/IKQWE+hEUHcGKv625DBwwcgDE2ljPjHnZ7n+DEiETmDJxD2+C26CQdbYPbMmfgHBIjEr27jGjSY2DNg33LPR2J15k+MBajTsNb33tptytBcCOpKrcOJUlqCWwBOqqqmnvD2CRgEkB0dHS3M2fOuC7KKlBVlf+8OpyQQ3dSEBDN+Nl3YPLTeyQWb5RmziEw1ERAsAnZqqDVSdW+Ha/IMgWbt5A2bVqZsaj58/Dr29flBU5kRcZqt2LUGkkvSMdX50uwT7BLr+EWSwZDXjrM/BU0DWftQWW8ud7M3I2prJ7Rm47NxIdloW6SJGmPqqrdq/I9lX7nlSTJH/gv8MSNyRdAVdWFqqp2V1W1e1hYWFViqDHZrmCxyiiqSs75I/il+JMbFMetg5uK5HsdW7Eda6HM2gUHa9Qd6RoP9AnWaXT46n25aLnI0K+G8lXqVy6/hlskPQZXzoC59vZN1xUT+95CkI+eN9aLWbDQsFQqAUuSpMeRfJerqupV73gWq8xG8yXuW7Cd2BfW8vGCKWh09+KjyaHLPQmeDs9rlNTEXvvewZp3R/qNJ/sER/pHkhCWwNpTdSShtb0LgqJhxwJPR+J1gnz0TO3fms3HMtlxsmE3sBAalsqsgpaAJcARVVXfdn9IlSfbFXacyOZwWg4fj+vI7r8l0Ti1G0U+oYTdEUf9rHFVPa7sjlTCk32CAYa0GsKxnGOcuFIH9pJqddBjEpzZ6ijOIZTycK+WhAca+fe6o25fUS8I3qIyM+A/AGOBgZIk7fvt1zA3x1UpVrvCLY00TIlMJeSz4ax6/jYsgXcQ4nOe+B7hFIuOK9e4ujsS4PE+wYNbDkYjaerOLLjLWOj/AvhHOPoFWwvA7mUVvDzEpNfy+G1t2Hv2Cj8ecd0iQUHwZpVZBb1VVVVJVdVOqqom/PZrTW0EdzMmrUqElInSPAF5/AYCLj2MImkZEvQWEdk78dGKT9Il3NUdyZN9gkN9Qrk1/FbWnV5XN2ZNWj00iYPlf4KXQ+GDwZC6wZGIBe7vHkXLEF/eWG/GrtSB/5+CUEN1uhpFkd3K9sKzjNv8JK/M+hMZgZ0Ia3UZw8QlmHa8C/ZiT4foNdzVHalUn+B9v9J8wXzH6udav4mKXwAAH3hJREFU6jA1tNVQzuSe4cjlI7VyvWqzy3ByM3w+DjIOgmJ3fF0xCk5tETNhQK/V8NQdcZgv5vHN/vOeDkcQ3K7OJmBZkcnIS6dTUDyfDVtJh6Z/Q6+3MTtyLsnFl5Bb90cy1L82g9VVbnekx2peE7ukT3Dm3Lmk3nY7qov3AFfk9ha3o5N0rDu1rtauWS32Ytj0qvOxja+KD4u/uTO+Ke2bBvL2hmNYxSMkoZ6rswkYq51mhmYEN2pKcXY+0Z0j6NZVy9JB81ly9FOs3ccj2VxTcrG+uL470pTfuiP5BhpcVhPbv28/sNsp+OUXl5yvMoKMQfRq1ou1p9eiqF78hq33gUspzscupTjGBTQaib8MiePc5UJW7jrr6XAEwa3qZAK2FcmcPXKFr/9zkAXTN/PtomOERPoSd2sUIUcz+EN4Eia/cNAaPR2q1ynpjoQE//33Hvb/lOayc/t07oQ2KIj8zVtcds7KGNJyCBkFGezP9OLVxbZCaFJOBdcm7R3jAgD92oSR2CqYd39MxWIVt+aF+qvOJWDFrpBmzmHt+4dL7WddvziFSxZ/pIzLjLzlXorsRY6tH4JTkiTRLK4xaeYcly1gkrRa/Hr3Jv/nn1GV2puNDoweiFFr9O7V0FojDHjB+diA58WHxetIksSzQ+KwWGX2n7tyrciOxSojV2O/uiB4qzqXgCvaz7prwwV8B9xGo6AmGMQb2k1FtW1MYa6VnHTXPbP179cXe1YWRSm1tyjKT+9H36i+rD+9Hlnx0hmTVget+sLolRDRCTQ6x9eRnzheFx8WS+nWIphV0/5Adr71WpGd+xZsZ5M5U8yKhXqjziXgm+1nNTYJQSks8u7i/F4iKq4x4KgP7Sp+vXuDJJG/eZPLzlkZQ1sN5XLRZXZl7KrV61aJwQ9iBsGj6+DFTEcyvpoGerFY8EayXeFUVgHTV/xKSnoudkUlJT2XCct280tqtpgJC/VCnUvAN93PailGMorZb2UEhvoQEGLivAsTsC44GFOnePK31O5z4D7N+jA9YTqxjWNRVAWLzeKds2GtzpGIJY1jD/C65+D8Hk9H5XWsdoV3fnBe5nT2D8dEkR2hXqhzCVirk0gcXs5+1mEt0Oo1bi8AUZ9ExTXm/LEcFBcWPvDv25eiAweRL1922TlvRlEVWgW1YsqGKXT9uCvj1o5j6/mtWGy1tyWqyjrcCzof0abQCZNei/lintMxc0YePnrRUUqo++pcAtZoNQQ0kRk8oUPZ/aztQtCaxOy3KprFNabYIpOdlu+yc/r37QeqSsHWrS47Z0VkRWZn+k6e2vwU5hwzdtWOOcfMjJ9mkJyR7J0zYQBTILQfAQf/K1ZB36DIZicuPMDpWFxEAIW2anbwEgQvUucSsGJXWPPiUrLP53H3E52ZMrc/9z7djRYdQ9CbxMy3qqLaNqb7sJb4Nzaiqiq2Ynu1OiNd7/+3d+fxVZV34sc/zznnLrkXkpCQELMQlgRCwp4YAwKirdUyrbi0U1dspxWVRUX7m9+oM/P7o/NzOi/7+vXXiYG61GmtFsvUvcLUFUHLIggCgcQga9gSCATIzV3Ovc/8cUPH6g255C7nXvK8Xy9fSe6TnOfr4eZ88+zOqkryHnwQ97RpyFCIkMdDyExcEvQH/TRsbYhY1rC1AX8w/scixs3kW8HXCU1vWh1JSrHrGg9ePSZi2ZKvj8FhpN2jS1G+Iu3exa88+zhntWp2/vE1nG4HQhPYHHpMOzkNZHanQW6Rm9d/sTU+ZwQTPiXJXjqcA3fNp2nCRPbdehtda9YQStAOWU7Dye5TuyOW7T65G6eRuGMRYzZiFmSVqG7oLzF0jelluTwzr4aqwkwMTVBVmMnT82qYXpaLoX7flYtAWr2LpZR0rjmLFBpfW3St1eGkvXNnBP/p6ca4nREcMk261q3n0P0P4GtqgmAQX1MTrQsW0rV+fUJawl7TS1l2WcSysiFleM3EHosYE02DSbfA5++HZ0Qrf+GyG8wem8eKu6fR/C/f5Je3V5OdYcNlVz1dysUhrRLw68/9Oz7n5WQEtzF8UpXV4aS9RJwRjN9Pe319xKL2+ifAH//uYLtuZ9GURRHLFk5eiF1P7LGIMZt8KyDh0xetjiTlGLqG22Gga4KHVnzK/S9uUXtEKxeNtEnAUkqOv9tGSDO44p6rrA7nopCIM4KF04mvJfLyEV9LC8IZ/+5gQzOoLail/qp6KnIqMIRBRU4F9VfVU1tQm/prwnNGQumMcDd0OhyraJF7Z4/mcKeX1z89bHUoihIXKf5k+h+rVjyN3zGDDHMbo2setDqci8K5NdXHI8yAPndG8IUe1CC9Xhzl5eHu5y9xlJcjvV5EAo4qdNlczCiaQW1BLQ7dwZGuI7gMF6502eRi8q3w2gI4uAGG11kdTUqaPTaPioLBPPnB59w4pQhN698RmoqSKtKmBXxk5V5Cmo3pP5hudSgXjYScEWy3k3ff4ohFeYsXgT1x3cGGZuCyuTh89jDffPmbrNqXwntDf1nlXLC5YcvzVkeSsoQQ3Dt7NC1tZ3m3qc3qcBQlZmmRgN959TkCtplkBLZTMUO1DuKl1zOCF/T/jGDNMHDX1VG8tAHHuHFgGDjGjaN4aQPuurqkbJJSklnCiMwRrG1dm/C64sYxCKquh8ZXwB95WEAJnxdcPCSDpat3x+0QEUWxSlp0Qe9/dSdB5ze49G+rrQ7lonPujOBLRmdjd+r4vUHsztiWdWkuF+5Zs3DX1SHsdszjx9GHDk3qDmUzi2fy+6bf4wl40qgb+rbwOPCuN2DSzVZHk5IMXWP+rFH882uNbNzbwWWjcq0OSVH6LeVbwB+s+gMBYyZO33YmXDXL6nAuSpquoemCZYtWs311a1zWVGuGgeZy0fb443x+zbWQxOMJAWYVz8If8rPx6Mak1huT0ukwZITqhu7Dd6tLyHXbWfbB51aHoigxSfkEvHvFBoJGBlO+W2F1KBc1u9NgyDAXbftOx/W6GdXVSJ8P747GuF63L9X51bgMF2tak3soREyECLeC962Fk/utjiZlZdh1fnD5CFY3t7PrSHzfr4qSTCmdgNe//yamNgundwdT51xjdTgXvfwRmRzbfyauY2uu6vCwgWfzprhdMxo23ca0wmmsPbQ2vcYKJ90CCPh0udWRpLQ76kbgtuv8UrWClTSW0gm48berMW1uxl83wupQBoT80ky6T/s5e9IXt2saubnYR46ke1Pyj9ybWTSTo11He92mMiVll8DIWbD1d0nvtk8nWS4bt142nDc+PcyBEyl84pWinEfKJuDN697DFDNxendx2Y3XWR3OgJA/Inz6TNv++HbruWqq8WzZgkxyQplZPBMgvbqhAabcDqf2w/6PrI4kpf1wxih0TfD02j1Wh6Io/ZKyCXjrr97EtA1i7DXDrA5lwBhaHF6K1LYv8jms/eWqqSF0+nSvO2QlSr4rn4qcCtYeSqPlSAAV34Ir/xGGjgEZCi9LCqbokYoWKshycuOUYlZsOsjxs/HrtVGUZEnJBLz9kz9jhmbi9DYz45bvWB3OgGHYdHKLB8W9BZxRXQOA5+PkjgNDuBt6a9tWOn2dSa+7/yTkjYEXboKfDIVnr4Hdb6v1wRHMv2IU/mCIX3+0z+pQFOWCpWQC3rDsJUx7JqOujHwgt5I4+aWDadt/BhmK38QlW1EhRkFB0idiQXg5UlAGWXdkXdLr7pegCXs+gBXz4Oh2CAXDH5ffDHvXqJbwl4zOG8S1VQU8t24fZ7wBq8NRlAvSZwIWQjwrhGgTQuxIZCBmMITHb9K2t5FQcAZObwtX3nl7IqtUIsgfkYm/26SzvTtu1xRC4KqupnvT5qTPSJ4wdAJZjqz02RUr6IPVj0Uue/+xcLnyV+65YjSnvSbLNx6wOhRFuSDRtIB/DST08F2P3+Skx4PEx+fLGwnYs6i6thiPX/21n2zDK3O44/9OIysvg4Av2K8zgSNx1VRjtrcTOHgwLteLlq7pXF54OR8e+pCQTINZxbYMaNsZuaxtZ7hc+SuTSrKZPjqXZ9buxWcGrQ5HUaLWZwKWUq4BOhIVgBkMYRAky+bAZXMx/p7ruPyGIszpgwlJL2acEoDSt4AvSPuBM6xatp1li1bz8uOb2d94goAv9odaxrn1wFYsRyqeSYe3g8bjyd0MpF8C3ZBfGbksvzJcrnzFvbNH03bGx6tbDlkdiqJEzfIx4JA/SGvjKV752RaWLVzNyqXbyMzLZnJONZ+2bSEYUuM6yRAKhmht6mDlsu0cbz2LDEmOt55l5dLttDZ3xNwSdpSVoWVlWTIOPKNwBgKRHrOhdQdc+WjksisfCZcrXzGjbCjjizJ58oM9BOM4f0FREiluCVgIMV8IsUkIsam9vT2qnwkGTFqbOlj1ZONfPfRXPdVIW8tZMm1ZBFHd0MkQNCUb39gbsWzjG3sJmrE91ISmMeyRR8hbtAgZChHyeAiZyfm3zXZm8091/8RN5TcRkiE8AQ9mKEXfV7oR3ojjlhehYCJoRvjjzcvDr+tpcX5K0gkhuPeKMvYc7+KtxqNWh6MoUYlbApZSPiWlrJFS1uTl5UX1M6EgfPxm5D1vP165nzFZFTgNNeaVDIZd48ThyMtcOg51Ydhie6uEPB60DCcH711A04SJ7Lv1NrrWrCHkSfwuRp6AhyHOISx8dyFTfzuVeavm8eGhD/EEUnQHJbsbyq6Gv/sv+Mc2+N7zkDEk/LrSq2vHFzAi18WyDz5Pr+1HlQHL0i5ow6Gf96Fvd9jwmmrMKxlMf4jcwsgP+JwiN2ag/13QIdOka916Dt13P76mJggG8TU10bpgIV3r1ye0JWyGTDYc2cCS1UtoPtlMUAZpPtnM4vcWs/HoxtRuCdvdoOnw4m3wzv+xOqKUp2uC+bNGs621kz9/fsLqcBSlT9EsQ1oOrAPGCiFahRA/jFflpi943od+wG9iF7Z4Vaech24Iaq8bGbGs9tsj0Q3R/4v7/bTX10csaq9/Avz+/l+7r6qDfhq2NkQsa9jagD+YuLrjpmouHNwApw9bHUnKu3FqEXmDHeqQBiUtRDML+hYp5SVSSpuUslhK+au4Va5D7d+URiyrnVOK0CSGYY9Xdcp5aLpG8dgc5iyYwNCS8JaUQ0sGMWfBBIrH5sR0RrBwOnvdhtLX0oJwOvt97b44DWevhzHsPrkbp5G4uuOm8vrwx52vWxtHGnDadH44YyRrW46zvTWddj9TBiJLu6B1m8Gw0Zlcc1fVXz/0766iqCIHm10l32SyOXRKq3K58cfV3P3EbObcO5HSqlxsDj2m60qvF0d5ecQyR3k50uuN6frn4zW9lGWXRSwrG1KG10xc3XEztBzyq2Dna1ZHkhZuu2w4g52GagUrKc/yZUgvPVzPiUNnuX7JFO55YjY3PDSV4eNzsWeormcraLqGzaGzaeU+nnv0zwSDcZjMYreTd9/iiEV5ixdBAv/Qsut2Fk1ZFLFs4eSF2PU0+SOvci4cWAdn1Azfvgx22vjB5SMYlefmjDdASEo8flPtKaCkHEsT8P7PdtLlq2LHSy/jcNkQmsDuNNBtaqmF1XIucYOEk0diPwBAMwzcdXUUL23AMW4cGAaOceMoXtqAu64OzUjcv7ehGdQW1FJ/VT0VORUYwqAip4L6q+qpLajF0NLkvVY5F5Cw6w2rI0kLP5oxirHDBvO9J9dT/ugqvrNsHaub29XuekpKsfTp894vnidofJ3iarXUKNXkFoUnx3Uc7iK/NDPm62kuF+5Zs3DX1SHsdoInT6INGZLQ5HuOy+ZiRtEMagtqcegOjnmOkZeRh01Po16W/ArIqwh3Q9feZXU0Kc0Mhtiw9wSLlm/5y2s7j5zmR89t4pl5Ncwem4cRw5wGRYkXy96Fvu5uTM9EHN17+da9kbsIFetk5WWgGYKOXpaJ9YdmGGguFwfnz+fg3fckJfmeY2gGLpuLlXtXcs1L17CnMw0Pca+cC/s/grNtVkeS0vzBED9/O/Kkv5+/8xk+U3VFK6nBsgT80k//Db8jH2fBPqtCUM5D0zWGDHP3uk47Fs6qKrwtLYQSuPyoNzXDwmcTf3z046TXHbPKuSBDqhu6D06bTvOxMxHLmo+eIcMW26RCRYkXyxJw9748bP6TzH34IatCUPqQU+im48jZuF/XWTUeAgF8zZ/F/dp9uWTQJRQPKk7PBJxfCbnlsPNVqyNJad5AkLHDIp8lPrZgMN0BdWKSkhosScDvPP8feDPGYdg+YXB2thUhKFHIKXRztsOHvzu+E1ec46sA8DZaczrRpQWXsunYpvQ4nvCLhAi3gvd9CF3HrY4mZdl1jQevHhOxbMnXx+Aw1PivkhoseSfuf+cwIhRg+vzrrKheidK5Xco64jAT+otsRUVoWVl4G3fE9brRurTgUk77T9NyMvI4YUqruj7cDd30R6sjSVmGrjG9LJdn5tVQVZiJoQmqCjN56o5qppflqglYSspI+izolq2b8Nuqcfg/oaL64WRXr1yAnML/mQldMCorbtcVQpBRVUW3RS3gL44Dj80Za0kM/TZsPOSMgsZXofr7VkeTslx2g9lj85g2OpcMm87hzm62HDiFy54my86UASHpfwqu/eXLhHQ7JbOjOzFJsU5mbgaGTYvrTOhznFVV+D5rIeTzxf3afUnrceBz3dB714Cnw+poUpqha7gdBpomePbDfSz5/VbaTqfBzmfKgJHUBNx1uhPTNxWn5zO+8f0fJbNqpR+EJsgpdHPicAImYo0fD6aJ77PkT8SCNB4HhvDe0DKouqEvwB3TSjFDkuUbD1odiqL8RVIT8Cs//RkBRw4Zw9V2euki5xJ3wlrAYO1ErLQdB75kEmSXqr2hL8DIoW5mjcnjdxv3E1BbUiopIqkJ2HeoBJvvODc8/L+SWa0Sg4rphdz0v6uRUhLwBQnF6eFlKypEz86me4c1E7HSej3wuW7oPauh+6TV0aSNeXWlHDvt463GY1aHoihAEhPwyqeX4c0ow8j4lAx35DOAldQS8AXxdflZtWw7yxau5uXHN7O/8QQBX+zrKIUQ4Q05GnfGIdILl9bjwBCeDR0yoWml1ZGkjSsr8ikeksFz6/ZZHYqiAElMwEc/7EQL+rhy0feSVaUSg1AwRGtTB6ue3MHx1rPIkOR461lWLt1Oa3NHXFrCzqoqfC3WTMSCNB8HLpwKWcNVN/QF0DXB7XWlbNjbQfPRyDtlKUoyJSUBb1+3Bq9jCo7AZkZWTkxGlUqMgqZk4xt7I5ZtfGMvQTP2Ywqd46vCE7Gam2O+Vn+k9TiwEFB5HXz+HnjVwfPR+l5NCQ5DU61gJSUkJQF//OxbSM3GqDkjk1GdEgeGXet1H+iOQ10YttjfOhk9E7HUOHA/VV4PoQA0r7I6krQxxG3n25MKeWXLIU57A1aHowxwCU/AnSeOEwhW4/TsZPbf3pbo6pQ4Mf2hv+yE9WU5RW7MQOzdtkZhIXlLlpB59dXIUIiQx0PITN55rWk/DlxUDVf/BEZdEd4dy98FQXXebV/mTSvF4w/y0uZWq0NRBriEJ+DX/vXnmPYsBpWr2ZrpRDcEtddF7rGo/fZIdEPEXIfs7sZeWsqBu+bTNGEi+269ja41awh5PDFfO1ppPQ5sdsOQUnjhu/CTofDsNbD77XAiVno1sTibySXZ/Hb9fqSMfShFUforoQnYDATwt4/G7jvGjX+vlh6lE03XKB6bw5wFExhaMghNEwwtGcScBRMoHpuDFuN+uiHTpGvdeg498AC+piYIBvE1NdG6YCFd69cnrSWctuPAQRP2fAAr5sHR7RAKhj8uvzm8S5ZqCZ/XvGml7Gnv4qPdJ6wORRnAErox6pvLnsCXMQm38Q42xy2JrEpJAJtDp7Qql8KybGwOHb83iN2px5x8AfD7aa+vj1jUXv8E7ro6MBK/b++lwy7l7ol3k+/KJyRDeE0vdt2OoaX4nsFBH6x+LHLZ+4/ByFmgp/j/g4XmTLiEf3lzF79Zt48Z5UOtDkcZoBLaAj6xOYhuevjGQz9IZDVKAmm6RtCULFu4mub1R+OTfAHhdOJridzq9LW0IJzOuNTTl0xHJuVDyrnrrbuY+tupzFs1jw8PfYgnkLxu8H6xZUBbL2uo23aGy5VeOW06N19awru7jnHoVLfV4SgDVMIS8Cfv/QmvYxK20GYKR5YnqholCTIG27A7dTrb4/egkl4vjvLI7wtHeTnSm/hN882QyYYjG/jxBz+m+WQzQRmk+WQzi99bzMajGzFDKdyNG+iG/MrIZfmV4XLlvG6rKwXghfX7LY5EGagSloC3/u4jpBBU3qDW/aY7IQRZ+S462+LYKrTbybtvccSivMWLwG6PX1298Af9NGxtiFjWsLUBf9Cf8Bj6TXfAlY9GLrvykXC5cl5F2Rl8bdwwXvz4IN5A7Lu7KcqFSkgCNgN+ArIGZ/cOpn37hkRUoSRZVl4Gp+LYAtYMA3ddHcVLG3CMGweGgWPcOIqXNuCuq0NLwviv03Cy+9TuiGW7T+7GaSSnG7xfdCM8znvLi1AwETQj/PGWF9X47wW4c9oIOrr8rNx+xOpQlAEoIb+lnUfbMW2DyB1jzRaDSvxl5Wfw+ZZ2gsEQepzGgTWXC/esWWRMnoyemUnI60U4nUlJvgBe00tZdhnNJ7+6E1fZkDK8pheXzZWUWPrF7oayq8MJ13BCZ2v4o13ttR6ty8tyGZXn5jfr9nPj1GKrw1EGmIS0gHNLirjzsenMfWBJIi6vWCArz4UMSc4cj+/YrGYY+D5roWn8BLzbtiUt+QLYdTuLpiyKWLZw8kLseuK7wWOmG+GEe/oQ/GIibF9hdURpRQjBHXWlfHrwFNtaT1kdjjLARJWAhRDXCiGahRC7hRD/0Nf3tx84w5tLt9G68yT+brXd28UgKz88qzaeE7HOcYwJT8bqbVZ0ohiaQW1BLfVX1VORU4EhDCpyKqi/qp7agtrUX4r0RdnDIW8ctLxldSRp56bqYlx2nefWqclYSnL1+YQRQuhAA3A10Ap8LIR4XUrZ6zlyUhI+OefJRubcXUVJVQ6G3Ra/qJWky84Pd8V2tnuA3Lhe28jJQc/NxZvkBAzgsrmYUTSDqflTcdlcdAW6cNvc6ZV8zxnzDVjXAN7T4My0Opq0kem0ccOUIv5zcyuPzBlHjjsNej6Ui0I0LeBaYLeUco+U0g+8CMyNtoKNK/cTCsW+baFirYzBNmwOnVNtiVne4igvT3oL+BxDMzA0g5rna3hu53PpmXwByq8JnxG8532rI0k786aNwG+GWLHpoNWhKAOI6GsvVCHEd4BrpZQ/6vn6DuAyKeWiL33ffGA+gNuRWZ0zuKDndcgbPpjNmzdvTkD86W4ocNzqINKAuk/RU/cqOuo+RU/dq+iMlVIOvpAfiNuf+lLKp4CnAIQQm856O2vide2LlRBik5RS3ac+qPsUPXWvoqPuU/TUvYqOEGLThf5MNF3Qh4CSL3xd3POaoiiKoij9FE0C/hgoF0KMFELYgZuB1xMblqIoiqJc3PrsgpZSmkKIRcCfAB14VkrZ2MePPRWP4AYAdZ+io+5T9NS9io66T9FT9yo6F3yf+pyEpSiKoihK/CX0OEJFURRFUSJTCVhRFEVRLBDXBHyhW1YOVEKIEiHE+0KInUKIRiHE/VbHlMqEELoQYosQ4o9Wx5KqhBDZQog/CCGahBC7hBDTrI4pVQkhlvT83u0QQiwXQqTwsVfJI4R4VgjRJoTY8YXXcoQQbwshWno+DrEyxlTRy716vOf3b5sQ4hUhRHZf14lbAv7ClpXfBCqBW4QQvZwYPuCZwENSykqgDlio7tV53Q/ssjqIFPcL4L+klBXAJNT9ikgIUQTcB9RIKccTnlh6s7VRpYxfA9d+6bV/AN6VUpYD7/Z8rUS+V28D46WUE4HPgIf7ukg8W8AxbVk5kEgpj0gpP+n5/Azhh2WRtVGlJiFEMfA3wDNWx5KqhBBZwCzgVwBSSr+UUh3t0zsDyBBCGIALOGxxPClBSrkG6PjSy3OB3/R8/hvg+qQGlaIi3Ssp5VtSSrPny/WE98w4r3gm4CLgixuptqKSSp+EECOAKcAGayNJWf8f+HsgZHUgKWwk0A78R09X/TNCCHUocARSykPAz4ADwBGgU0qpjpDq3TAp5ZGez48Cw6wMJo38HbCqr29Sk7AsJIQYBLwEPCClPG11PKlGCPEtoE1KqfYRPz8DmAosk1JOAbpQXYUR9YxhziX8R0sh4BZC3G5tVOlBhtesqnWrfRBCPEp4mPGFvr43nglYbVl5AYQQNsLJ9wUp5ctWx5OiLgeuE0LsIzykcZUQ4nlrQ0pJrUCrlPJcL8ofCCdk5au+DuyVUrZLKQPAy8B0i2NKZceEEJcA9HxsszielCaE+D7wLeA2GcUmG/FMwGrLyigJIQTh8bpdUsr/Z3U8qUpK+bCUslhKOYLw++k9KaVqrXyJlPIocFAIMbbnpa8BvZ7XPcAdAOqEEK6e38OvoSasnc/rwJ09n98JvGZhLClNCHEt4eGy66SUnmh+Jm4JuGfw+dyWlbuAFVFsWTlQXQ7cQbhFt7XnvzlWB6WktcXAC0KIbcBk4DGL40lJPb0EfwA+AbYTfgaqrRYBIcRyYB0wVgjRKoT4IfBT4GohRAvh3oOfWhljqujlXj0BDAbe7nmm/7LP66itKBVFURQl+dQkLEVRFEWxgErAiqIoimIBlYAVRVEUxQIqASuKoiiKBVQCVhRFURQLqASsKIqiKBZQCVhRFEVRLPDf4AU3DDDXNhEAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 576x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots(1, 1, figsize=(8, 4))\n",
"\n",
"# We want to evaluate the system on 30 linearly\n",
"# spaced times between t=0 and t=3.\n",
"t = np.linspace(0., 3., 30)\n",
"\n",
"# We simulate the system for different values of k.\n",
"for k in np.linspace(0., 1., 5):\n",
" # We simulate the system and evaluate $v$ on the\n",
" # given times.\n",
" v = spi.odeint(func1, v0, t, args=(k,))\n",
" # We plot the particle's trajectory.\n",
" ax.plot(v[:, 0], v[:, 1], 'o-', mew=1, ms=8,\n",
" mec='w', label=f'k={k:.1f}')\n",
"ax.legend()\n",
"ax.set_xlim(0, 12)\n",
"ax.set_ylim(0, 6)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"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.7.7"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment