HW_2_V.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import scipy.integrate as integ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "#rk-4 time stepping\n",
    "#inputs-present solution:yp,present time:tp,f-derivative function,time-step-dt,present time:tp\n",
    "#returns next timestep tn and next solution yn\n",
    "\n",
    "def rk4(yp,tp,f,dt):\n",
    "    \n",
    "    #rk-4 coefficients\n",
    "    \n",
    "    k1=f(yp,tp)\n",
    "    \n",
    "    k2=f(yp+dt*(k1/2),tp+(dt/2))\n",
    "    \n",
    "    k3=f(yp+dt*(k2/2),tp+(dt/2))\n",
    "    \n",
    "    k4=f(yp+dt*k3,tp+dt)\n",
    "    \n",
    "    #next time\n",
    "    \n",
    "    tn=tp+dt \n",
    "    \n",
    "    #next solution\n",
    "    \n",
    "    yn=yp+(dt/6)*(k1+2*k2+2*k3+k4)\n",
    "    \n",
    "    return tn,yn\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "#a rk-4 solver\n",
    "\n",
    "#inputs-initial solution y0,step size dt,final time tf,function f\n",
    "\n",
    "def rk4_solver(y0,dt,tf,f):\n",
    "    \n",
    "    #number of derivatives\n",
    "\n",
    "    N=len(y0)\n",
    "\n",
    "    #number of time steps\n",
    "    \n",
    "    t=np.arange(0,tf,dt)\n",
    "\n",
    "    Nt=len(t)\n",
    "    \n",
    "    #solution initialization\n",
    "    \n",
    "    y=np.zeros((Nt,N))\n",
    "\n",
    "    y[0]=y0\n",
    "\n",
    "    #time counter\n",
    "\n",
    "    tp=0 \n",
    "\n",
    "    #rk4 time stepping\n",
    "\n",
    "    for i in range(1,Nt):\n",
    "\n",
    "        #solution at time-step i\n",
    "\n",
    "        y[i]=rk4(y[i-1],tp,f,dt)[1]\n",
    "\n",
    "        #updating time\n",
    "\n",
    "        tp=rk4(y[i-1],tp,f,dt)[0]\n",
    "\n",
    "    return y\n",
    "\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "#numerical solution\n",
    "\n",
    "tf=2*np.pi\n",
    "\n",
    "#rhs function\n",
    "\n",
    "def f(y0,tp):\n",
    "    return np.array([y0[1],-y0[0]])\n",
    "\n",
    "#initial y0\n",
    "\n",
    "y0=[1,0]\n",
    "\n",
    "def numerical_solution_error(dt):\n",
    "\n",
    "    numerical_solution=rk4_solver(y0,dt,tf,f)[0:,:1] #first column of solution\n",
    "\n",
    "    #number of time points\n",
    "\n",
    "    t=np.arange(0,tf,dt)\n",
    "\n",
    "    Nt=len(t)\n",
    "\n",
    "    #reshaped numerical solution\n",
    "\n",
    "    numerical_solution=numerical_solution.reshape(Nt)\n",
    "\n",
    "    #exact solution\n",
    "\n",
    "    exact_solution=np.cos(t)\n",
    "\n",
    "    #error \n",
    "\n",
    "    error=np.max(np.abs(numerical_solution-exact_solution))\n",
    "\n",
    "    return error\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "#numerical error for standard library function\n",
    "\n",
    "def library_error(dt):\n",
    "\n",
    "    #number of time points\n",
    "\n",
    "    t=np.arange(0,tf,dt)\n",
    "\n",
    "    Nt=len(t)\n",
    "\n",
    "\n",
    "    numerical_solution=integ.odeint(f,y0,t)[0:,:1] #first column of solution\n",
    "\n",
    "    \n",
    "    #reshaped numerical solution\n",
    "\n",
    "    numerical_solution=numerical_solution.reshape(Nt)\n",
    "\n",
    "    #exact solution\n",
    "\n",
    "    exact_solution=np.cos(t)\n",
    "\n",
    "    #error \n",
    "\n",
    "    error=np.max(np.abs(numerical_solution-exact_solution))\n",
    "\n",
    "    return error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.18554080283684e-14\n"
     ]
    }
   ],
   "source": [
    "#checking my code accuracy\n",
    "\n",
    "dt=0.001\n",
    "\n",
    "print(numerical_solution_error(dt))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.183776390944473e-08\n"
     ]
    }
   ],
   "source": [
    "#checking library code accuracy\n",
    "\n",
    "dt=0.001\n",
    "\n",
    "print(library_error(dt))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkgAAAHJCAYAAAB+GsZPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABz7ElEQVR4nO3dd3gU9drG8e+mF0joJZBQpJOEEoqgVOm9WI4eFWyv2BCRIyAdpCiCyFHsApZzjg2QpghKb0qJBBJKIAjSa0ICabvz/jEQTYA0EmaT3J/ryqW/2ZnZZ5cs+/CbuWdshmEYiIiIiEgaF6sLEBEREXE2apBEREREMlCDJCIiIpKBGiQRERGRDNQgiYiIiGSgBklEREQkAzVIIiIiIhmoQRIRERHJQA2SiIiISAZqkERE8tHhw4ex2WzMmzcvx9tGRkYyfvx4Dh8+nOd1iUjm1CCJiDipyMhIJkyYoAZJxAJqkESkwEpJSSE1NfWGj12+fPmW9m0YBleuXLmlfYhIwaUGSUQsd+DAAR566CHKlSuHp6cndevW5d133023zpo1a7DZbHz++ee8/PLLVKpUCU9PT6Kjoxk4cCDFihUjIiKCTp06Ubx4ce655x4Azp8/z7PPPkulSpXw8PCgevXqjBo1iqSkpHT7t9lsPP/887z//vvUrVsXT09P5s+ff9Oaq1atSo8ePVi4cCGhoaF4eXlRvXp1Zs+ena3XvGHDBu655x6KFy+Oj48PLVu2ZNmyZWmPz5s3j/vuuw+Adu3aYbPZcn2oTkRyzs3qAkSkaIuMjKRly5YEBQUxY8YMKlSowIoVKxg8eDBnz55l3Lhx6dYfOXIkLVq04P3338fFxYVy5coBkJycTK9evXj66acZMWIEqampJCYm0q5dOw4ePMiECRMIDQ1l/fr1TJ06lfDw8HQNCcCiRYtYv349Y8eOpUKFCmn7vpnw8HCGDBnC+PHjqVChAl9++SUvvvgiycnJDBs27KbbrV27lo4dOxIaGsonn3yCp6cnc+bMoWfPnvz3v//lgQceoHv37kyZMoVXX32Vd999l8aNGwNwxx135OZtFpGcMkRELNS5c2ejcuXKRmxsbLrlzz//vOHl5WWcP3/eMAzDWL16tQEYrVu3vm4fAwYMMADj008/Tbf8/fffNwDj66+/Trf89ddfNwDjp59+SlsGGP7+/mnPl5UqVaoYNpvNCA8PT7e8Y8eOhp+fn5GQkGAYhmHExMQYgDF37ty0de68806jXLlyxqVLl9KWpaamGsHBwUblypUNh8NhGIZhfPPNNwZgrF69Ols1iUje0SE2EbFMYmIiP//8M3379sXHx4fU1NS0n27dupGYmMiWLVvSbdO/f/+b7i/jY7/88gu+vr7ce++96ZYPHDgQgJ9//jnd8vbt21OyZMls11+/fn0aNGiQbtlDDz1EXFwcO3bsuOE2CQkJbN26lXvvvZdixYqlLXd1deWRRx7hzz//ZN++fdmuQUTyhxokEbHMuXPnSE1N5d///jfu7u7pfrp16wbA2bNn021TsWLFG+7Lx8cHPz+/6/ZfoUIFbDZbuuXlypXDzc2Nc+fOZWvfN1OhQoWbLsu472suXLiAYRg3fK6AgIBMtxWR20fnIImIZUqWLJk2c/Lcc8/dcJ1q1aqlG2dsdjJbXrp0abZu3YphGOkeP336NKmpqZQpUyZb+76ZkydP3nRZ6dKlb7hNyZIlcXFx4cSJE9c9dvz4cYDr6hKR208zSCJiGR8fH9q1a8fOnTsJDQ2lSZMm1/3crNHIjnvuuYf4+HgWLVqUbvlnn32W9vit2LNnD7///nu6Zf/5z38oXrx42knVGfn6+tK8eXMWLFiQ7jICDoeDL774gsqVK1OrVi0APD09AXS5ARELaAZJRCz19ttvc/fdd9OqVSueeeYZqlatyqVLl4iOjmbJkiX88ssvud73o48+yrvvvsuAAQM4fPgwISEhbNiwgSlTptCtWzc6dOhwS7UHBATQq1cvxo8fT8WKFfniiy9YuXIlr7/+Oj4+PjfdburUqXTs2JF27doxbNgwPDw8mDNnDrt37+a///1v2kxWcHAwAB9++CHFixfHy8uLatWq3VLTKCLZowZJRCxVr149duzYwaRJkxg9ejSnT5+mRIkS1KxZM+08pNzy8vJi9erVjBo1iunTp3PmzBkqVarEsGHDrrt8QG40bNiQxx57jHHjxnHgwAECAgKYOXMmL730UqbbtWnThl9++YVx48YxcOBAHA4HDRo0YPHixfTo0SNtvWrVqjFr1izefvtt2rZti91uZ+7cuWknmYtI/rEZhmFYXYSISEFTtWpVgoODWbp0qdWliEg+0DlIIiIiIhmoQRIRERHJQIfYRERERDLQDJKIiIhIBmqQRERERDJQgyQiIiKSga6DlEsOh4Pjx49TvHjxHN+eQERERKxhGAaXLl0iICAAF5ebzxOpQcql48ePExgYaHUZIiIikgtHjx6lcuXKN31cDVIuFS9eHDDf4Ix3EBcRERHnFBcXR2BgYNr3+M2oQcqla4fV/Pz81CCJiIgUMFmdHqOTtEVEREQyUIMkIiIikoEOseUzu91OSkqK1WVIIeXh4ZFpCkNERHJHDVI+MQyDkydPcvHiRatLkULMxcWFatWq4eHhYXUpIiKFihqkfHKtOSpXrhw+Pj66VpLkuWvX4jpx4gRBQUH6HRMRyUNqkPKB3W5Pa45Kly5tdTlSiJUtW5bjx4+TmpqKu7u71eWIiBQaOnkhH1w758jHx8fiSqSwu3ZozW63W1yJiEjhogYpH+mQh+Q3/Y6JiOSPIt8gXb58mSpVqjBs2DCrSxEREREnUeQbpMmTJ9O8eXOryxAREREnUqQbpAMHDrB37166detmdSmShYEDB9KnTx+ryxARkSLCaRukdevW0bNnTwICArDZbCxatOi6debMmUO1atXw8vIiLCyM9evX5+g5hg0bxtSpU/OoYhk/fjwNGza0ugwRESkM7KmWPr3TNkgJCQk0aNCAd95554aPf/XVVwwZMoRRo0axc+dOWrVqRdeuXTly5EjaOmFhYQQHB1/3c/z4cb7//ntq1apFrVq1slVPUlIScXFx6X6kYEtOTr5umd1ux+Fw5Hhfud1OREQyMAz4/X8wuyGcj7GyDucHGAsXLky3rFmzZsagQYPSLatTp44xYsSIbO1zxIgRRuXKlY0qVaoYpUuXNvz8/IwJEybcdP1x48YZwHU/sbGx16175coVIzIy0rhy5YphGIbhcDiMhKQUS34cDke23g/DMIwffvjBuOuuuwx/f3+jVKlSRvfu3Y3o6Oh06xw9etR44IEHjJIlSxo+Pj5GWFiYsWXLFmPu3LnXvTdz5841YmJiDMDYuXNn2j4uXLhgAMbq1asNwzCM1NRU4/HHHzeqVq1qeHl5GbVq1TJmzZqV7nkHDBhg9O7dO9P6N27caLRq1crw8vIyKleubLzwwgtGfHx82uNVqlQxJk2aZAwYMMDw8/MzHn30UWPu3LmGv7+/sWTJEqNu3bqGq6urcejQIeP8+fPGI488YpQoUcLw9vY2unTpYuzfvz9tXzfb7nbL+LsmIlKgnTtkGJ/1MYxxfubPkpfy/CliY2Nv+v39dwXyQpHJycls376dESNGpFveqVMnNm3alK19TJ06Ne3w2rx589i9ezdjx4696fojR45k6NChaeO4uDgCAwOz9VxXUuzUG7siW+vmtciJnfHxyN4fc0JCAkOHDiUkJISEhATGjh1L3759CQ8Px8XFhfj4eNq0aUOlSpVYvHgxFSpUYMeOHTgcDh544AF2797Njz/+yKpVqwDw9/fn1KlTWT6vw+GgcuXKfP3115QpU4ZNmzbxf//3f1SsWJH7778/W7VHRETQuXNnJk2axCeffMKZM2d4/vnnef7555k7d27aetOnT2fMmDGMHj0agA0bNnD58mWmTp3Kxx9/TOnSpSlXrhwPPfQQBw4cYPHixfj5+TF8+HC6detGZGRk2gUZb7SdiIjkgj0VtsyB1VMg9Qq4ekKbV+CuFy0rqUA2SGfPnsVut1O+fPl0y8uXL8/Jkyfz5Tk9PT3x9PTMl307i/79+6cbf/LJJ5QrV47IyEiCg4P5z3/+w5kzZ/jtt98oVaoUADVq1Ehbv1ixYri5uVGhQoUcPa+7uzsTJkxIG1erVo1Nmzbx9ddfZ7tBmj59Og899BBDhgwBoGbNmsyePZs2bdrw3nvv4eXlBUD79u3TXdJhw4YNpKSkMGfOHBo0aACQ1hht3LiRli1bAvDll18SGBjIokWLuO+++wCu205ERHLh+E5YPBhO7jLHVVtBz7eh9B2WllUgG6RrMl4kzzCMXF04b+DAgXlU0Y15u7sSObFzvj5HZs+dXQcPHmTMmDFs2bKFs2fPpp1Tc+TIEYKDgwkPD6dRo0ZpzVFeev/99/n444/5448/uHLlCsnJyTk64Xv79u1ER0fz5Zdfpi0zDAOHw0FMTAx169YFoEmTJtdt6+HhQWhoaNo4KioKNze3dJd/KF26NLVr1yYqKuqm24mISA4kJ5gzRlvmgOEArxLQ6TVo9DA4wUVwC2SDVKZMGVxdXa+bLTp9+vR1s0rOwGazZfswl5V69uxJYGAgH330EQEBATgcDoKDg9NOZvb29s7xPl1czByAYRhpy67diuWar7/+mpdeeokZM2bQokULihcvzvTp09m6dWu2n8fhcPD0008zePDg6x4LCgpK+39fX9/rHvf29k7XWP+91r/L2IBn3E5ERLIpehUsfQkuXg1W1e8HXV+HYs5zqoLzf2vfgIeHB2FhYaxcuZK+ffumLV+5ciW9e/e2sLKC69y5c0RFRfHBBx/QqlUrwDz89HehoaF8/PHHnD9//oazSB4eHtfdE6xs2bIAnDhxgkaNGgEQHh6ebp3169fTsmVLnn322bRlBw8ezFH9jRs3Zs+ePekO+eVWvXr1SE1NZevWrWmH2M6dO8f+/fvTZqJERCQXEs7CjyMh4mtz7B8I3WdALWuOsmTGaWP+8fHxhIeHp32ZxsTEEB4enhbjHzp0KB9//DGffvopUVFRvPTSSxw5coRBgwZZWHXBVbJkSUqXLs2HH35IdHQ0v/zyS7qT0gEefPBBKlSoQJ8+fdi4cSOHDh3iu+++Y/PmzQBUrVo17c/p7NmzJCUl4e3tzZ133sm0adOIjIxk3bp1aSdIX1OjRg22bdvGihUr2L9/P2PGjOG3337LUf3Dhw9n8+bNPPfcc4SHh6edR/TCCy/k+L2oWbMmvXv35qmnnmLDhg38/vvvPPzww1SqVEkNuIhIbhgGhP8H3mliNkc2F7jzWXh2i1M2R+DEDdK2bdto1KhR2qzD0KFDadSoUVrS7IEHHmDWrFlMnDiRhg0bsm7dOpYvX06VKlWsLLvAcnFx4X//+x/bt28nODiYl156ienTp6dbx8PDg59++oly5crRrVs3QkJCmDZtGq6u5nlO/fv3p0uXLrRr146yZcvy3//+F4BPP/2UlJQUmjRpwosvvshrr72Wbr+DBg2iX79+PPDAAzRv3pxz586lm03KjtDQUNauXcuBAwdo1aoVjRo1YsyYMVSsWDFX78fcuXMJCwujR48etGjRAsMwWL58eVqCTUREsun8Ifi8Dyx6Bq5cgPLB8OQq6DIVPItZXd1N2YybnXAhmYqLi8Pf35/Y2Fj8/PzSPZaYmEhMTEzaVb5F8ot+10TEadlTYPO7sGaaGd1384I2w6HlC+Bq3T82M/v+/rsCeQ6SiIiIOLFjO8zo/qkIc1ytNfSYZXl0PyfUIImIiEjeSIo3o/tb3zOj+94lodNkaPiQU0T3c0INkoiIiNy6Ayth6VCIvRrdD74XukyDYmWtrSuX1CCJiIhI7sWfgR9HwO5vzbF/EPSYCTU7WlvXLVKDJCIiIjlnGBD+JawYBYkX/4rutx3p1Om07FKDJCIiIjlz7iAsHQIx68xxhRDoORsqNba0rLykBklERESyx54Cm/4Na1+H1EQzut92JLR4ztLofn5QgyQiIiJZ+3M7LBkMp3ab4+ptocdbUKq6pWXlFzVIIiIicnNJl+CXybD1fcAA71LQeQo0+EeBi+7nhNPeakQkv61ZswabzcbFixdztX3VqlWZNWtWntYkIuJU9q+Ad+80r2uEASH3w/O/QcMHC3VzBGqQpIC51aZGRESyIf40fPMY/Od+iPsTSgTBw99B/4/At4zV1d0WOsQmIiIiJsOAnZ/DT6MhMfav6H67V8HD1+rqbivNIN0OhgHJCdb85OBexIZh8MYbb1C9enW8vb1p0KAB3377bdpjHTp0oEuXLly7v/HFixcJCgpi1KhRwF+zO8uWLaNBgwZ4eXnRvHlzIiIi0j3Ppk2baN26Nd7e3gQGBjJ48GASEhLSHk9KSuKVV14hMDAQT09PatasySeffMLhw4dp164dACVLlsRmszFw4MAsa79m+fLl1KpVC29vb9q1a8fhw4ezfE/Gjx9PUFAQnp6eBAQEMHjw4Juue+TIEXr37k2xYsXw8/Pj/vvv59SpU+n21bBhQz744AMCAwPx8fHhvvvuu242bO7cudStWxcvLy/q1KnDnDlzsqxTROSWnTsI83vC4hfM5qhCKDy1GjpPLnLNEWgG6fZIuQxTAqx57lePZ/sXe/To0SxYsID33nuPmjVrsm7dOh5++GHKli1LmzZtmD9/PiEhIcyePZsXX3yRQYMGUb58ecaPH59uP//61794++23qVChAq+++iq9evVi//79uLu7ExERQefOnZk0aRKffPIJZ86c4fnnn+f5559n7ty5ADz66KNs3ryZ2bNn06BBA2JiYjh79iyBgYF899139O/fn3379uHn54e3t3e2aj969Cj9+vVj0KBBPPPMM2zbto2XX3450/fj22+/5a233uJ///sf9evX5+TJk/z+++83XNcwDPr06YOvry9r164lNTWVZ599lgceeIA1a9akrRcdHc3XX3/NkiVLiIuL44knnuC5557jyy+/BOCjjz5i3LhxvPPOOzRq1IidO3fy1FNP4evry4ABA7L15ygikiOpybBpNqx9A+xJ4OZtzhjd+Sy4Ft02oei+ckknISGBmTNn8ssvv9CiRQsAqlevzoYNG/jggw9o06YNlSpV4oMPPuCRRx7h1KlTLFmyhJ07d+Lunv7aF+PGjaNjR/MS8/Pnz6dy5cosXLiQ+++/n+nTp/PQQw8xZMgQAGrWrMns2bNp06YN7733HkeOHOHrr79m5cqVdOjQIa2Oa0qVKgVAuXLlKFGiRLZrf++996hevTpvvfUWNpuN2rVrExERweuvv37T9+TIkSNUqFCBDh064O7uTlBQEM2aNbvhuqtWrWLXrl3ExMQQGBgIwOeff079+vX57bffaNq0KQCJiYlp7wnAv//9b7p3786MGTOoUKECkyZNYsaMGfTr1w+AatWqERkZyQcffKAGSUTy3p/bzBmj05HmuHq7q9H9atbW5QTUIN0O7j7mTI5Vz50NkZGRJCYmpjU21yQnJ9OoUaO08X333cfChQuZOnUq7733HrVq1bpuX9eaFDAbmtq1axMVFQXA9u3biY6OTpsxAXP2xeFwEBMTQ0REBK6urrRp0ybbLzE7tUdFRXHnnXdi+1vq4u913sh9993HrFmzqF69Ol26dKFbt2707NkTN7frPzZRUVEEBgamNUcA9erVo0SJEkRFRaU1SEFBQWnN0bUaHA4H+/btw9XVlaNHj/LEE0/w1FNPpa2TmpqKv79/tt8PEZEsJV2CnyfBrx8CBviUhs5TIfT+Qp9Oyy41SLeDzeb0x28dDgcAy5Yto1KlSuke8/T0TPv/y5cvs337dlxdXTlw4EC293+tMXE4HDz99NM3PJcnKCiI6OjofKndyMG5WNcEBgayb98+Vq5cyapVq3j22WeZPn06a9euvW7WzDCMdM1XVsuvufaYzWZLex0fffQRzZs3T7eeq6trjusXEbmhfT/Aspch7pg5bvAgdJoMvqWtrcvJqEESwJzt8PT05MiRI5nO3rz88su4uLjwww8/0K1bN7p370779u3TrbNlyxaCgoIAuHDhAvv376dOnToANG7cmD179lCjRo0b7j8kJASHw8HatWvTDrH9nYeHBwB2uz1HtderV49FixZdV2dWvL296dWrF7169eK5556jTp06RERE0Lhx+vsN1atXjyNHjnD06NG0WaTIyEhiY2OpW7du2npHjhzh+PHjBASY56Rt3rwZFxcXatWqRfny5alUqRKHDh3in//8Z5a1iYjkyKVT8MMrELnIHJeoAj1nwR3tM9uqyFKDJAAUL16cYcOG8dJLL+FwOLj77ruJi4tj06ZNFCtWjAEDBrBs2TI+/fRTNm/eTOPGjRkxYgQDBgxg165dlCxZMm1fEydOpHTp0pQvX55Ro0ZRpkwZ+vTpA8Dw4cO58847ee6559JOPo6KimLlypX8+9//pmrVqgwYMIDHH3887STtP/74g9OnT3P//fdTpUoVbDYbS5cupVu3bnh7e2er9kGDBjFjxgyGDh3K008/zfbt25k3b16m78m8efOw2+00b94cHx8fPv/8c7y9valSpcp163bo0IHQ0FD++c9/MmvWrLSTtNu0aUOTJk3S1vPy8mLAgAG8+eabxMXFMXjwYO6//34qVKgAmEm3wYMH4+fnR9euXUlKSmLbtm1cuHCBoUOH3voftIgUPQ4H7PwMVo69Gt13Ne+d1nYkeGTvNIwiyZBciY2NNQAjNjb2useuXLliREZGGleuXLGgstxzOBzG22+/bdSuXdtwd3c3ypYta3Tu3NlYu3atcfr0aaN8+fLGlClT0tZPSUkxmjVrZtx///2GYRjG6tWrDcBYsmSJUb9+fcPDw8No2rSpER4enu55fv31V6Njx45GsWLFDF9fXyM0NNSYPHly2uNXrlwxXnrpJaNixYqGh4eHUaNGDePTTz9Ne3zixIlGhQoVDJvNZgwYMCDL2q9ZsmSJUaNGDcPT09No1aqV8emnnxqAceHChRu+HwsXLjSaN29u+Pn5Gb6+vsadd95prFq1Ku3xKlWqGG+99Vba+I8//jB69epl+Pr6GsWLFzfuu+8+4+TJk2mPjxs3zmjQoIExZ84cIyAgwPDy8jL69etnnD9/Pt3zfvnll0bDhg0NDw8Po2TJkkbr1q2NBQsW3LDGgvq7JiK3yZn9hvFpV8MY52f+vN/aMI6HZ71dIZbZ9/ff2QwjFydnCHFxcfj7+xMbG4ufn1+6xxITE4mJiaFatWp4eXlZVOHtt2bNGtq1a8eFCxfSEmbyl/Hjx7No0SLCw8PzbJ9F9XdNRLKQmgwbZ8G66WBPNgM77UZB80FFOroPmX9//13RfpdEREQKm6O/wuLBcMZMD1OjA3SfCSWvPz1Abk4NkoiISGGQGAc/T4TfPsaM7peBLtMg5F5F93NBh9hySYfYxBnod01EANi7DJYNg0tXr7nX8J/Q6TXwKWVtXU5Ih9hEREQKu0snYfm/IGqxOS5ZzYzuV29rZVWFghokERGRgsbhgB3zYeU4SLoa3W/5ArQdAe7eVldXKKhBEhERKUjO7IclL8KRTeY4oDH0mg0VQqytq5BRgyQiIlIQpCbDhrdg/ZtXo/u+0H40NH8aXHQ7orymBklERMTZHdlizhqd2WuOa3aC7jOgRJC1dRViapBEREScVWIsrJoA2z4xx75lzeh+cH9F9/OZi9UFiPNo27YtQ4YMSRtXrVqVWbNmpY1tNtt1N3wVEZF8ErUE3m3+V3PU6GF47ldd1+g20QyS3NRvv/2Gr6+v1WWIiBQtcSdg+TDYu9Qcl6oOPd+Gaq2trauIUYMkN1W2bNlb3kdKSgru7u453i45ORkPD49bfv7cPpfdbsdms+HikrNJ1txuJyKCwwHb58Kq8ZAUBy5u0HIwtHlF0X0L6G/x28AwDC6nXLbk51YulJ7xEBvAiRMn6Nq1K97e3lSrVo1vvvkm7bHDhw9js9n4+uuvadu2LV5eXnzxxRecO3eOBx98kMqVK+Pj40NISAj//e9/0+23bdu2PP/88wwdOpQyZcrQsWNHHn/8cXr06JFuvdTUVCpUqMCnn35607o3bdpE69at8fb2JjAwkMGDB5OQkJDudb322msMHDgQf39/nnrqKebNm0eJEiVYunQp9erVw9PTkz/++IMLFy7w6KOPUrJkSXx8fOjatSsHDhxI29fNthMRyZHTe2FuV1g21GyOKoXB/62FDuPUHFlEM0i3wZXUKzT/T3NLnnvrQ1vxcffJs/2NGTOGadOm8fbbb/P555/z4IMPEhwcTN26ddPWGT58ODNmzGDu3Ll4enqSmJhIWFgYw4cPx8/Pj2XLlvHII49QvXp1mjf/632ZP38+zzzzDBs3bsQwDM6fP0/r1q05ceIEFStWBGD58uXEx8dz//3337C+iIgIOnfuzKRJk/jkk084c+YMzz//PM8//zxz585NW2/69OmMGTOG0aNHA7BhwwYuX77M1KlT+fjjjyldujTlypXjoYce4sCBAyxevBg/Pz+GDx9Ot27diIyMTJsZu9F2IiLZkpoE62fC+hngSDGj+/eMhWZPKbpvMTVIkiP33XcfTz75JACTJk1i5cqV/Pvf/2bOnDlp6wwZMoR+/fql227YsGFp///CCy/w448/8s0336RrkGrUqMEbb7yRbrvatWvz+eef88orrwAwd+5c7rvvPooVK3bD+qZPn85DDz2UdrJ5zZo1mT17Nm3atOG9995Lu19Z+/bt09W0YcMGUlJSmDNnDg0aNABIa4w2btxIy5YtAfjyyy8JDAxk0aJF3HfffQDXbSciki1/bIYlg+HsfnNcqwt0exNKBFpblwBqkG4Lbzdvtj601bLnzkstWrS4bhweHp5uWZMmTdKN7XY706ZN46uvvuLYsWMkJSWRlJR03QngGbcDePLJJ/nwww955ZVXOH36NMuWLePnn3++aX3bt28nOjqaL7/8Mm2ZYRg4HA5iYmLSZrpu9FweHh6EhoamjaOionBzc0vXxJUuXZratWsTFRV10+1ERDJ15aJ5ntH2q7PavuWg6+tQv6/SaU5EDdJtYLPZ8vQwl7OxZfhAZ2x8ZsyYwVtvvcWsWbMICQnB19eXIUOGkJycnOl2AI8++igjRoxg8+bNbN68mapVq9KqVaub1uJwOHj66acZPHjwdY8FBf11QbUbPZe3t3e613Kz87cMw0i3XsbtRERuyDDMm8oufwXiT5rLGj8KHSeCd0lra5PrqEGSHNmyZQuPPvpounGjRo0y3Wb9+vX07t2bhx9+GDCbmAMHDqQ7b+lmSpcuTZ8+fZg7dy6bN2/msccey3T9xo0bs2fPHmrUqJGNV5O5evXqkZqaytatW9MOsZ07d479+/dnq3YRkTRxx2HZMNi3zByXrgE9ZkG1m/+DT6ylBkly5JtvvqFJkybcfffdfPnll/z666988sknmW5To0YNvvvuOzZt2kTJkiWZOXMmJ0+ezHaT8eSTT9KjRw/sdjsDBgzIdN3hw4dz55138txzz/HUU0/h6+tLVFRU2rlSOVGzZk169+7NU089xQcffEDx4sUZMWIElSpVonfv3jnal4gUUQ6HeaHHVRMg+ZIZ3b/7JWg1DNy9rK7OqR09f5nAUtYdfVHMX3JkwoQJ/O9//yM0NJT58+fz5ZdfUq9evUy3GTNmDI0bN6Zz5860bduWChUq0KdPn2w/Z4cOHahYsSKdO3cmICAg03VDQ0NZu3YtBw4coFWrVjRq1IgxY8akpeByau7cuYSFhdGjRw9atGiBYRgsX748V9d2EpEi5nQUfNrZvOhj8iWo3BSeXm/eYFbN0U1dSkxhzKLdtJm+mt8On7esDptxKxfKKcLi4uLw9/cnNjYWPz+/dI8lJiYSExNDtWrV0lJTknuXL18mICCATz/99Lp0XFGn3zURJ5SSaMb2N7xlRvc9isE946DpE4ruZ+GnPScZ+/0eTsYlAjCkQ02GdKiVp8+R2ff33+kQmzgth8PByZMnmTFjBv7+/vTq1cvqkkREMnd4Iyx5Ec5dvaBs7W7QbTr4V7a2Lid3Oi6RcYv38MNu8+T1KqV9mNo3hJY1ylhWkxokcVpHjhyhWrVqVK5cmXnz5uHmpl9XEXFSVy7CyrGwY745LlYeur4B9Xorup8Jh8Pgf78dZeoPUVxKTMXVxcb/ta7Oi/fUxMvd2tk2feOI06pateot3SpFRCTfGQZELoIfhkP8KXNZ4wHQcYKi+1k4eCaekQsi+DXGPM+oQWV/pvYLpV7AzQ973U5qkERERHIj9k8zur//B3Ncuib0fBuq3mVtXU4uOdXB+2sP8s4v0STbHfh4uPJyp9oMbFkVVxfnmW1Tg5SPHA6H1SVIIacZNhELOOzw2yfw8wRIjgcXd2g1FO4eqnRaFrb/cYGRC3ax/1Q8AG1rl+W1PsFULul8F1NWg5QPPDw8cHFx4fjx45QtWxYPDw9daVnynGEYnDlzBpvNpssOiNwup/aYJ2H/+Zs5rtwMes2Gcrp4bGbik1KZ/uNePtvyB4YBpX09GNuzHr0aBDjt96MapHzg4uJCtWrVOHHiBMePH7e6HCnEbDYblStXxtVV0WGRfJWSCOumw8ZZ4EgFj+LQYRw0eQJcdEnBzKyKPMWY73dzItaM7t8bVplR3epS0tfD4soypwYpn3h4eBAUFERqaip2u93qcqSQcnd3V3Mkkt9i1puzRucPmuPa3a9G9ytZW5eTO30pkQmLI1kWcQIwo/tT+oZwl4XR/ZxQg5SPrh360OEPEZEC6MoF+GkM7PzcHBerYDZG9XRNtsw4HAZfbTvK1OVRxF2N7j/Vyozue3sUnH/QqUESERH5O8OAPQvghxGQcNpcFvYYdBgP3iWsrMzpZYzuh1b2Z2q/EOoH+FtcWc6pQRIREbnm4lFY9jIcWGGOy9SCnrOhSgtr63JyyakOPlh7kH+vjiY51YG3uysvd6rFwJZVcXMtmOdoqUESERFx2OHXD+HnSZCSYEb3Ww+Du18CN0+rq3NqO45cYOR3Eew7dQmANrXM6H5gKeeL7ueEGiQRESnaTu6GJYPh2HZzHHinGd0vW9vaupxcQYzu54QaJBERKZpSrsDaN2DTbDO67+lnnmcU9pii+1nIGN3v37gyo7s7f3Q/J9QgiYhI0XNoLSwdAucPmeO6Pc2by/oFWFqWs8sY3Q8qZUb3765ZMKL7OaEGSUREio7L583ofvgX5rh4Rej2JtTtYW1dTs4wDL767ShT/hbdf7JVNYbcU6tARfdzokg3SDExMTz++OOcOnUKV1dXtmzZgq+vr9VliYhIXjMM2P0d/DgCEs4ANmj6BNwzFrwKXgT9djp0Nbq/9Wp0P6SSGd0PrlS437ci3SANHDiQ1157jVatWnH+/Hk8PZVUEBEpdC4egaVDIXqlOS5bx4zuBzW3ti4nl5zq4MN1B5n9S+GJ7udEkW2Q9uzZg7u7O61atQKgVKlSFlckIiJ5ymGHrR/AL6+Z0X1XD2g1DO4eouh+FnYeucCIv0X3W9cqy+RCEN3PCadtAdetW0fPnj0JCDDjgosWLbpunTlz5lCtWjW8vLwICwtj/fr12d7/gQMHKFasGL169aJx48ZMmTIlD6sXERFLnYyAj++BFSPN5iioJQzaCG2HqznKRHxSKuMX76Hfe5vYd+oSpXw9mPVAQ+Y/1rRINUfgxDNICQkJNGjQgMcee4z+/ftf9/hXX33FkCFDmDNnDnfddRcffPABXbt2JTIykqCgIADCwsJISkq6btuffvqJlJQU1q9fT3h4OOXKlaNLly40bdqUjh073rCepKSkdPuKi4vLo1cqIiJ5JuUKrJkGm/4Nhh08/aHjBGg8QNH9LPwcdYoxi3Zz/Gp0v1/jSozuXo9ShSi6nxM2wzAMq4vIis1mY+HChfTp0ydtWfPmzWncuDHvvfde2rK6devSp08fpk6dmuU+N2/ezIQJE/jxxx8BmD59OgD/+te/brj++PHjmTBhwnXLY2Nj8fPzy8nLERGR/HBwNSx9CS7EmON6vc3ofvEK1tbl5E5fSmTCkkiW7TKj+4GlvJnSN4RWNctaXFn+iIuLw9/fP8vv7wLZTicnJ7N9+3Y6deqUbnmnTp3YtGlTtvbRtGlTTp06xYULF3A4HKxbt466devedP2RI0cSGxub9nP06NFbeg0iIpJHLp+Hhc/A533M5qh4APzjv3D/Z2qOMmFG94/QYcZalu06gauLjadbV+enIW0KbXOUE057iC0zZ8+exW63U758+XTLy5cvz8mTJ7O1Dzc3N6ZMmULr1q0xDINOnTrRo8fNr4Ph6emplJuIiDMxDIj41ozuXz4L2KDZU9B+DHhpZj8zMWcTGLlgF1sOmdH94Ep+TOsXWuij+zlRIBukazLe68UwjBzd/6Vr16507do1r8sSEZH8duEPWDYUoleZ47J1zfunBTazti4nl2J38OG6Q7z984G06P7QjrV47K6iEd3PiQLZIJUpUwZXV9frZotOnz593aySiIgUIvZU2PoerJ4CKZfN6H7rV+CuF8GtaJ5MnF07j1xg5III9p40o/utapZhSt+QIpdOy64C2SB5eHgQFhbGypUr6du3b9rylStX0rt3bwsrExGRfHPid1g8GE6Em+Mqd0PPWVCmppVVOb34pFTeXLGP+ZsPYxhQyteDMT3q0qdhpRwddSlqnLZBio+PJzo6Om0cExNDeHg4pUqVIigoiKFDh/LII4/QpEkTWrRowYcffsiRI0cYNGiQhVWLiEieS74Ma6bC5nfN6L6XP3ScBI0eUXQ/C7/sPcXohX+L7jeqxOgeRTe6nxNO2yBt27aNdu3apY2HDh0KwIABA5g3bx4PPPAA586dY+LEiZw4cYLg4GCWL19OlSpVrCpZRETy2sFfYMkQuPiHOa7fF7q8DsV1OkVmzlxKYsKSPSz9W3R/cp8QWtdSOi27CsR1kJxRdq+jICIiuZBwDla8Crv+Z479KkH3GVBbwZrMGIbBN9v+ZPLyKGKvpOBigydbVWdIh5r4eDjtnMhtld3vb71bIiLiPAwDdn1t3iLk8jnABs2fhvajwbO41dU5tZizCby6IILNh84BUD/Aj9f7K7qfW2qQRETEOZyPMaP7B38xx+Xqm9H9yk2srcvJXYvuz/75AEmpDrzcXRjasRaP31VN0f1boAZJRESsZU+FLXPM6H7qFXD1hDZXo/uu7lZX59TCj15kxHe70kX3J/cJIai0ovu3Sg2SiIhY53g4LH4BTu4yx1VbQc+3ofQdlpbl7BKSUnnzp33M22RG90v6uDOmRz36NlJ0P6+oQRIRkdsvOcGcMdoyBwwHeJWATq9Bo4dBX/CZWr33NKMX7ebYxSsA9G1UidHd61K6mG6HlZfUIImIyO0VvQqWvgQXj5jj4P7QZRoUK2dtXU7uzKUkJi6NZMnvxwGoXNKbyX1DaKPofr5QgyQiIrdHwln4cSREfG2O/QPN6H6tztbW5eQMw+Cb7X8yedlf0f0n7q7GSx1rKbqfj/TOiohI/jIM+P2/5nWNrlwAmws0HwTtRoFnMaurc2qHzybw6sIINh38K7o/rV8oIZUV3c9vapBERCT/nD9kHk47tMYclw82o/uVwiwty9ml2B18tP4Qb6/6K7r/UodaPHG3ovu3ixokERHJe/YU895pa6aZ0X03L2gzHFq+oOh+Fn4/epERCyKIOhEHwN01yjClr6L7t5saJBERyVvHdsDiwXAqwhxXaw09Zim6n4WEpFRm/LSfeZticFyN7o/uXo9+jRXdt4IaJBERyRtJ8WZ0f+t7ZnTfuyR0mgwNH1J0Pwur951m9MK/ovt9GgYwpkc9RfctpAZJRERu3YGVsHQoxF6N7ofcB52nQjFF0DNzNj6JiUsiWazovtNRgyQiIrkXfwZ+HAG7vzXH/kHQYybU7GhtXU7OMAy+3f4nk5dHcfGyGd1//K5qDO2k6L6z0J+CiIjknGFA+JewYhQkXjSj+3c+C21HKrqfhcNnExi1KIKN0WZ0v15FP6b1DyG0cglrC5N01CCJiEjOnDsIS4dAzDpzXCEEes6GSo0tLcvZpdgdfLw+hlmr9pOU6sDTzYWXOprRfXdF952OGiQREckeewps+jesfR1SE8HNG9qNNGeOFN3P1K4/LzL8u7+i+3fVKM2UviFUKe1rcWVyM2qQREQka39uhyWD4dRuc1y9LfR4C0pVt7QsZ5eQlMrMlfuZu9GM7pe4Gt3vr+i+01ODJCIiN5cUD7+8BlvfBwzwLgWdp0CDfyi6n4U1+04z6m/R/d5Xo/tlFN0vENQgiYjIje3/CZYNhdij5jj0AbM58i1jbV1O7mx8EpOWRvJ9uBndr1TCm9f6BtOudjmLK5OcUIMkIiLpxZ+GH4bDngXmuESQeTitRgdr63IShmHcdPl3O47x2vJIYq9G9x9r+Vd0/2bb3S4G1j5/btiwWXYo0mZY/SdWQMXFxeHv709sbCx+fn5WlyMimF9QV1KvEJ8ST3xKPAnJCcSnxHM55fJfy1IS0j12bXxtnYSUBOKT40lxpFj9cnKsIH4BimTm1eav8mCdB/N0n9n9/tYMkohYLsWR8lcTk3yTpiXDY9etk5xAQmoCDsNh9csRkUJADZKI5IphGCTaE9NmXNKamBvMxqSbubm67HLq5bTtEu2JeVqbi80FX3dfirkX++u/Hr74uvlSzKPYDR+7Nr724+laRE6kTU2BbR/DpnfBngRuXnD3EGg8EFxds70bG855wnZ+Hp6JOh7HxKWR7D91CYBm1UoxqltdKpfyybwmZ32vnLAub3dvy55bDZJIEWN32ElITbjuEFPaLExKwnWzNDc6XJWQkoDdsOdpbZ6unumblwzNzM0eK+ZeDB93n7THvN28FaHOjj+3weIX4HSkOb6jvXmuUcmqlpbl7C4npzLzp/18mhbdL8mobnW5N6yyfu8KETVIIgWAYRgkO5KvP8SUfIMZm789lrHRSUhJ4ErqlTytzYYtffPikWF2JovHinkUw9fNXMddFxu8PZIuwc+T4NcPAQN8Sps3lg29X9H9LKzdf4ZRCyP484L5OerVIICxPRXdL4zUIInkI4fhuL55+dvMzXUNTPL1Dc21/091pOZpbR4uHn81Lzc57JTVIali7sXwcvPCxabbJBQY+36EZS9D3J/muMGD0Gky+Ja2ti4nd+5qdH+RovtFhhokkRtIsaekTz3dZDbmpo9dOxyVejnPa/v7eTI3nI3JxiEpX3dfPFw98rw2cWKXTsEPr0DkInNcsqp5OO2O9lZW5fTSovvLIrl4Nbo/sGU1Xu5UC19PfYUWZvrTlULjZhHv6w473eD8m4yNTl5HvN1sbjdvXnJwSMrH3UezNZIzDgfs/AxWjoXEWLC5QovnoO1I8Mj8ZOKi7o9zCYxauJsN0WcBqFOhOK/3D6VBYAlrC5PbQg2SWO5mEe+cHpJKSEnI8+vAeLt55zoJ9feZGw8XD528Kbff2QOw5EX4Y6M5rtgQes2Gig0sLcvZpdodfLwhhlmr9pOY4sDTzYUXO9TkqVbVcXfVP1CKCjVIkivXZmv+HtXO7JBUZicPJ9mT8rQ2V5tr1k1LVoekPHzxcfPBzUUfESmAUpNh49uwbroZ3Xf3gXajoPkgcNXvdGYi/oxlxIJd7DkeB0DLO0ozpW8IVcv4WlyZ3G76pBQxqY5UElISrr8AXyZXFr7Z4aq8viCfl6vXTU8Yvhbhzs7MjZerl2ZrpOg6+issHgxnosxxjQ7QfSaUrGJtXU7ucnIqb63czycbzOi+v7c7o7srul+UqUEqAAzDIMmelKsL8OV3xNvF5mJGtG8W387BzI1ma0RuQWIc/DwRfvsYM7pfBrpMg5B7Fd3Pwrr9Z3hV0X3JQN9ITuazPZ/x85GfrzvPJtXI24j3tQvyZTYbk27m5iYnE+uCfCJOYO9yM7p/yYyg0/Cf0Ok18CllbV1O7lx8Eq8ti2LhzmPA1eh+n2Da1VF0X9QgOZ3jCcfZcXrHDR+77oJ8ub2GjS7IJ1I4XDoJy/8FUYvNcclq0HMWVG9rZVVOzzAMFu48xqSlkVxQdF9uQr8JTqbXHb1oUr5J+pmbq02Qt5u3It4iYkb3d8yHleMg6Wp0/67B0GY4WHjvqoLgyLnLjFoUwfoDiu5L5tQgOZl6petRr3Q9q8sQEWd1Zr8Z3T+yyRwHNDaj+xVCrK3LyaXaHXy6MYaZKxXdl+xRgyQiUhCkJsOGt2D9m2BPBndfaD8amj8NLq5WV+fUdh8zo/u7j5nR/RbVSzOlXwjVFN2XTKhBEhFxdke2mLNGZ/aa45qdoPsMKBFkbV1O7nJyKrNWHeCTDTHYHQb+3u6M6l6X+xTdl2xQgyQi4qwSY2HVBNj2iTn2LWtG94P7K7qfhXX7zzBqUQRHz5vR/Z4NAhjbox5liyu6L9mjBklExBlFLTETapdOmONGD0PHSYruZ+F8QjKTlkamRfcD/L14rW8w7euUt7gyKWjUIImIOJO4E7B8GOxdao5LVYeeb0O11tbW5eQyRvdtNhjQoirDOtemmKL7kgv6rRERcQYOB2yfC6vGQ1IcuLjBXS9C638pup+Fo+cv8+rC9NH9qf1CaBRU0uLKpCBTgyQiYrXTe82TsI9uMceVwqDnbKgQbG1dTi5jdN/DzYUX76nJ/7VWdF9unRokERGrpCbB+pmwfgY4Uszo/j1jodlTiu5nQdF9yW9qkERErPDHZlgyGM7uN8e1ukC3N6FEoLV1ObkryXbeWrU/fXS/W13ua6LovuQtNUgiIrfTlYvmeUbb55pj33LQ9XWo31fR/SysP3CGVxf+Fd3vEVqRcT3rK7ov+UINkojI7RK52Izux580x40fhY4TwVsnE2fmfEIyry2LZMEORffl9lGDJCKS3+KOm43Rteh+6RrQYxZUa2VpWc7OMAy+Dz/OxKWRnE9IVnRfbiv9homI5BeHw7wK9qoJkHzJjO7f/RK0GgbuXlZX59SOnr/MqEW7Wbf/DKDovtx+apBERPLD6ShYPBj+/NUcV25qRvfL17O2LieXancwd+NhZq7cz5UUu6L7Yhk1SCIieSkl0Yztb3jLjO57FIN7xkHTJxTdz8LuY7GMXBBBxLFYAO6sXoopfUOoXraYxZVJUaQGSUQkrxzeaF7w8dwBc1y7G3SbDv6Vra3LyV1JtjNr1X4+vhrd9/NyY1T3utzfJFDRfbGMGiQRkVt15SKsHAs75pvjYuWh6xtQr7ei+1nYcOAsry6M4Mj5ywB0D63IuJ71KFdc52iJtdQgiYjklmFA5CL4YTjEnzKXhQ2EDhPAu4SFhTm/CwnJTPpbdL+ivxeTegfToZ6i++Ic1CCJiORG7J+wbBjs/8Ecl64JPd+GqndZW5eTU3RfCgr9NoqI5ITDDr99Aj9PgOR4cHGHVkPh7qGK7mfh6PnLjF60m7VXo/u1yxdnav8QGiu6L06oSDdIb731Fh9//DGGYdChQwfefvttnRAoIjd3ao95Evafv5njwObmrFG5utbW5eRS7Q7mbTrMjJ/+iu4Pbl+D/2t9Bx5uiu6LcyqyDdKZM2d455132LNnD+7u7rRu3ZotW7bQokULq0sTEWeTkgjrpsPGWeBIBY/i0GEcNHkCXPQFn5k9x2MZ8d1f0f3m1UoxtZ+i++L8imyDBJCamkpiYiIAKSkplCtXzuKKRMTpHN5wNbofbY5rd78a3a9kbV1O7kqynVk/7+fj9X9F91/tZkb3XVw0Uy/Oz2n/6bNu3Tp69uxJQEAANpuNRYsWXbfOnDlzqFatGl5eXoSFhbF+/fps779s2bIMGzaMoKAgAgIC6NChA3fccUcevgIRKdCuXIDvn4d53c3mqFgFuP9zePA/ao6ysDH6LF3eXscHaw9hdxh0D6nIqpfb8I9mQWqOpMBw2hmkhIQEGjRowGOPPUb//v2ve/yrr75iyJAhzJkzh7vuuosPPviArl27EhkZSVBQEABhYWEkJSVdt+1PP/2Et7c3S5cu5fDhw3h7e9O1a1fWrVtH69at8/21iYgTMwzYs9CM7iecNpc1eRw6jAcvf0tLc3YXEpJ5bVkU3+34E1B0Xwo2m2EYhtVFZMVms7Fw4UL69OmTtqx58+Y0btyY9957L21Z3bp16dOnD1OnTs1yn9988w1r1qzh3XffBWD69OkYhsErr7xyw/WTkpLSNVtxcXEEBgYSGxuLn59fLl+ZiDiVi0dh2ctwYIU5LlPLvH9aFZ2bmBnDMFj8+3EmLonk3NXo/qN3VmFY59oU93K3ujyRdOLi4vD398/y+9tpD7FlJjk5me3bt9OpU6d0yzt16sSmTZuytY/AwEA2bdpEYmIidrudNWvWULt27ZuuP3XqVPz9/dN+AgMDb+k1iIgTcdhhy/sw506zOXJxhzYjYNAGNUdZOHr+MgPn/saL/wvnXEIytcoX49tBLZnQO1jNkRRoOT7ElpqaipeXF+Hh4QQHB+dHTVk6e/Ysdrud8uXTT9uWL1+ekydPZmsfd955J926daNRo0a4uLhwzz330KtXr5uuP3LkSIYOHZo2vjaDJCIF3MndsGQwHNtujgPvhF6zoezN/8EkYHcYzN0Y81d039WFF9rX4Ok2iu5L4ZDjBsnNzY0qVapgt9vzo54cyXjNIsMwcnQdo8mTJzN58uRsrevp6Ymnp2eO6hMRJ5ZyBda+AZtmm9F9Tz/oOAEaD1R0PwuRx+MYsWAXu/40o/vNrkb371B0XwqRXJ2kPXr0aEaOHMkXX3xBqVKl8rqmLJUpUwZXV9frZotOnz593aySiMh1Dq2FpUPg/CFzXLcndJ0OfhUtLcvZJabYmbXqAB+tN9Npxa9G9x9QdF8KoVw1SLNnzyY6OpqAgACqVKmCr69vusd37NiRJ8XdjIeHB2FhYaxcuZK+ffumLV+5ciW9e/fO1+cWkQLs8nn4aQyEf2GOi1eEbm9C3R7W1lUAbIw+y6sLI/jj3GUAuoVUYHzP+pTz0+1VpHDKVYP09zRZfomPjyc6OjptHBMTQ3h4OKVKlSIoKIihQ4fyyCOP0KRJE1q0aMGHH37IkSNHGDRoUL7XJiIFjGHA7u/gxxGQcAawQdMn4J6xiu5n4UJCMpOXR/HtdjO6X8HPi0l9gumo6L4Uck4b81+zZg3t2rW7bvmAAQOYN28eYF4o8o033uDEiRMEBwfz1ltv3bbrGGU3JigiFrt4BJYOheiV5rhsHTO6H9Tc2rqc3I2i+4/cWYV/KbovBVx2v79vqUHavn07UVFR2Gw26tWrR6NGjXK7qwJHDZKIk3PYYesH8MtrkJIArh7Q+l9w1xBw87C6Oqf254XLjF60mzX7zgBQq3wxpvYLJaxKSYsrE7l12f3+ztUhttOnT/OPf/yDNWvWUKJECQzDIDY2lnbt2vG///2PsmXL5rpwEZFbdjICFr8Ax3ea46CW0PNtKFvL2rqcnN1hMG/TYWb8tI/LyWZ0//n2NRik6L4UQbn6jX/hhReIi4tjz549nD9/ngsXLrB7927i4uIYPHhwXtcoIpI9KVdg5Tj4oI3ZHHn6m43RwGVqjrIQeTyOfnM2MmlpJJeT7TSrWorlL7Zi8D011RxJkZSrQ2z+/v6sWrWKpk2bplv+66+/0qlTJy5evJhX9TktHWITcTKH1sCSIXAhxhzX6w1d34DiFaysyuklpth5++cDfLTuEKlXo/sju9blH00V3ZfCKV8PsTkcDtzdrz9Jz93dHYfDkZtdiojkzuXzsGIU/P4fc1w8ALrPgDrdrK2rANh0Nbp/+Gp0v2twBSb0UnRfBHLZILVv354XX3yR//73vwQEBABw7NgxXnrpJe655548LVBE5IYMAyK+NaP7l88CNmj2FLQfA16a1c3MxcvJTF4WxTd/i+5P7F2fTvU12yZyTa4apHfeeYfevXtTtWpVAgMDsdlsHDlyhJCQEL744ou8rlFEJL0Lf8CyoRC9yhyXrWvePy2wmbV1OTnDMFiy6wQTl+zhbLwZ3X+4eRVe6aLovkhGuWqQAgMD2bFjBytXrmTv3r0YhkG9evXo0KFDXtcnIvIXeypsfR9WT4aUy+DqCW3+BS1fVHQ/C8cuXmH0wghWX43u1yxXjGn9QwircvtvFyVSEOS4QUpNTcXLy4vw8HA6duxIx44d86MuEZH0TvwOiwfDiXBzXOVuM6FWpoalZTk7u8Ng/qbDvPm36P5z7WowqG11PN1crS5PxGnluEFyc3OjSpUq2O32/KhHRCS95MuwZipsfhcMu3lrkI6ToNEj4KL4eWaiTsQxYkEEvx+9CEDTqiWZ2i+EGuWKW1uYSAGQq0Nso0ePZuTIkXzxxReUKqXpWRHJJwd/MaP7F/8wx/X7QpfXobjuA5aZxBQ7s38+wIfXovuebozoVocHmwYpui+STblqkGbPnk10dDQBAQFUqVIFX1/fdI/v2LEjT4oTkSIq4RyseBV2/c8c+1WC7jOhdhdr6yoANh08y6sL/orud6lfgQm961Ne0X2RHMlVg9SnT588LkNEBDO6v+trWDESLp8DbND8aWg/Gjx1WCgzFy8nM2V5FF9vM6P75f08mdg7mM6K7ovkSq5O0gZ4/PHHCQwMzPOCRKSIOh9jRvcP/mKOy9U3o/uVm1hbl5MzDIOlu04w4Wp0H+DhO4N4pUsd/BTdF8m1XN1qpHjx4kRERFC1atV8KKlg0K1GRPKIPRW2zIHVUyD1ihndbzscWg4GV33BZ+bYxSuMWbSbX/aeBqBGuWJM6xdCk6o6N1TkZvL1ViP33HMPa9asYeDAgbmtT0QEjofD4hfg5C5zXLWVGd0vfYelZTk7u8Pgs82HeXPFPhKuRvefbXcHz7S9Q9F9kTySqwapa9eujBw5kt27dxMWFnbdSdq9evXKk+JEpJBKTjBnjLbMAcMBXiWg02vQ6GGwKWWVmb0n4xj+3V/R/SZVSjKtv6L7InktV4fYXDK59ojNZisS10jSITaRXIpeBUtfgotHzHFwf+gyDYqVs7YuJ5eYYuffvxzgg7V/RfeHd63DQ80U3RfJiXw9xOZwOHJdmIgUUQln4ceREPG1OfYPNKP7tTpZW1cBsPngOV5dGEHM2QQAOtcvz4RewVTwV3RfJL/kqkH6u8TERLy89CEVkZswDPj9v+Z1ja5cAJsLNB8E7UaBZzGrq3NqsZdTmLI8iq+2HQWgXHEzut8lWNF9kfyWqwbJbrczZcoU3n//fU6dOsX+/fupXr06Y8aMoWrVqjzxxBN5XaeIFETnD5mH0w6tMcflg83ofqUwS8tydoZhsCziBOMXR3I2PgmAfzYPYnhXRfdFbpdc3cho8uTJzJs3jzfeeAMPj7/uoB0SEsLHH3+cZ8WJSAFlT4ENs2BOS7M5cvOCDuPh/9aoOcrC8YtXeHL+Np7/z07OxidxR1lfvhnUgsl9Q9QcidxGuZpB+uyzz/jwww+55557GDRoUNry0NBQ9u7dm2fFiUgBdGwHLB4MpyLMcbU20OMtRfezYHcYfL75MNOvRvfdXW0827YGz7ZTdF/ECrlqkI4dO0aNGjWuW+5wOEhJSbnlokSkAEqKN6P7W98zo/veJaHTZGj4kKL7Wdh7Mo4R30UQ/rfo/tR+IdQsr+i+iFVy1SDVr1+f9evXU6VKlXTLv/nmGxo1apQnhYlIAXJgJSwdCrFXo/sh90HnqVCsrLV1ObnEFDvv/BLN+2sPKrov4mRy1SCNGzeORx55hGPHjuFwOFiwYAH79u3js88+Y+nSpXldo4g4q/gz8OMI2P2tOfYPgh4zoWZHa+sqALYcOserCyI4pOi+iFPKVYPUs2dPvvrqK6ZMmYLNZmPs2LE0btyYJUuW0LGj/mIUKfQMA8K/hBWjIPGiGd2/81lo9yp4+Ga5eVEWezmFqT9E8b/f/h7dr0+X4IoWVyYif5erK2mLrqQtRdi5g7B0CMSsM8cVQqDnbKjU2NKynJ1hGCyPOMm4xXvSovsPNQ9ieJc6+HsrnSZyu+TrlbRFpAiyp8Cmf8Pa1yE1Edy8od1Ic+bIVV/wmTl+8QpjFu3m572nAbijrC9T+4XSrFopiysTkZvJVYNUsmRJbDdIpdhsNry8vKhRowYDBw7kscceu+UCRcQJ/LkdlgyGU7vNcfW2ZnS/VHVLy3J2N4ruP9O2Bs8pui/i9HLVII0dO5bJkyfTtWtXmjVrhmEY/Pbbb/z4448899xzxMTE8Mwzz5CamspTTz2V1zWLyO2SFA+/vAZb3wcM8C4FXaZC6AOK7mdh38lLjFiwi51HLgIQdjW6X0vRfZECIVcN0oYNG3jttdfSXSQS4IMPPuCnn37iu+++IzQ0lNmzZ6tBEimo9v8Ey4ZCrHkyMaEPQOcp4FvG2rqcXMbofjFPN4Z3qc0/m1dRdF+kAMnVSdrFihUjPDz8uotFRkdH07BhQ+Lj4zl48CChoaEkJCTkWbHORCdpS6EVfxp+GA57FpjjEkHm4bQaHaytqwDIGN3vWK88E3vXp6K/t8WVicg1+XqSdqlSpViyZAkvvfRSuuVLliyhVCnzpMOEhASKF9dUskiBYRiw83P4aTQkxprR/RbPQduRiu5nIfZKCtN+iOK/v5qzbWWLezKxV326BFe44fmaIuL8ctUgjRkzhmeeeYbVq1fTrFkzbDYbv/76K8uXL+f9998HYOXKlbRp0yZPixWRfHLuICx5EQ6vN8cVQqHXvyGgoaVlOTvDMPhhtxndP3PJjO4/2CyIEV0V3Rcp6HJ9HaSNGzfyzjvvsG/fPgzDoE6dOrzwwgu0bNkyr2t0SjrEJoVCajJsmg1r3wB7khndbz8Kmj8DrroKSGZOxF5hzKI9rIo6BUD1sr5M7RtC8+qlLa5MRDKT3e9vXSgyl9QgSYH35zZY/AKcjjTHd7Q3zzUqWdXSspydw2HwxdY/eOPHfcQnpaZF959tewde7oruizi7fL9QpN1uZ9GiRURFRWGz2ahXrx69evXC1VV/QYg4taRL8PMk+PVDwACf0tBlmnmDWZ0vk6n9py4x4rtd7Lga3W8cVIJp/UMV3RcphHLVIEVHR9OtWzeOHTtG7dq1MQyD/fv3ExgYyLJly7jjjjvyuk4RyQv7foBlL0PcMXPc4EHoNBl8dVgoM4kpduasjua9tQdJsSu6L1IU5OoQW7du3TAMgy+//DIttXbu3DkefvhhXFxcWLZsWZ4X6mx0iE0KlEun4IdXIHKROS5Z1Tycdkd7K6sqELYeOsfIhREcOqPovkhhkK+H2NauXcuWLVvSmiOA0qVLM23aNO66667c7FJE8oPDATs/g5Vjr0b3XaHl89BmBHj4WF2dU1N0X6Roy1WD5OnpyaVLl65bHh8fj4eHxy0XJSJ54OwBM7r/x0ZzXLEh9JoNFRtYWpazMwyDH3efZKyi+yJFWq4apB49evB///d/fPLJJzRr1gyArVu3MmjQIHr16pWnBYpIDqUmw8ZZsG462JPB3Qfaj4ZmTyu6n4XjF68w9vu/RffL+DK1n6L7IkVRrv62nD17NgMGDKBFixa4u5v/okpJSaF3797MmjUrL+sTkZw4+issHgxnosxxjQ7QfSaUrGJtXU7O7jD4fPNhpq/YR0KyHXdXG4Pa3MFz7Wooui9SROWqQSpRogTff/890dHRREVFYRgG9erVu+7ebCJymyTGwc8T4bePMaP7ZaDr6xDcX9H9LESdiGPEggh+P3oRgLAqJZnaL0TRfZEiLtsN0tChQzN9fM2aNWn/P3PmzFwXJCI5tHcZLBsGl46b44b/hE6vgU+pzLcr4hJT7Lz98wE+WneIVIdBcU83hnetw0PNghTdF5HsN0g7d+7M1npKd4jcJpdOwvJ/QdRic1yyGvScBdXbWllVgbDhwFlGLYrgj3OXAegaXIHxvepT3s/L4spExFlku0FavXp1ftYhItnlcMCO+bByHCRdje7fNRjaDAd3XZsnM+cTknltWSQLdpgXyqzg58XE3vXpVL+CxZWJiLNRpEWkIDmz34zuH9lkjgMam9H9CiHW1uXkDMNg4c5jTFoayYXLKdhsMKBFVV7uVIviXorui8j11CCJFASpSbDhLVg/42p03xfuGQPN/g9clLLKzB/nEhi1cDcbos8CUKdCcab2C6FRUEmLKxMRZ6YGScTZHdliRvfP7jPHNTtB9xlQIsjaupxcit3Bx+tjmLVqP0mpDjzdXHixQ02ealUdd1cXq8sTESenBknEWSXGwqoJsO0Tc+xb1ozu1++n6H4Wwo9eZMR3u9h70rzi/101SjO5TwhVy/haXJmIFBRqkEScUdQSM6F26YQ5bvQwdJyk6H4W4pNSeXPFPuZvPoxhQEkfd0Z3r0e/xpWUsBWRHFGDJOJM4k7A8mGwd6k5LnWHGd2v1trSsgqCVZGnGPP9bk7EJgLQr1ElRnWvS+linhZXJiIFkRokEWfgcMD2ubBqPCTFgYsb3PUitP6XovtZOB2XyPgle1gecRKAoFI+TO4bTKuaZS2uTEQKMjVIIlY7vdeM7h/dYo4rhUHP2VAh2Nq6nJzDYfDf344w7Ye9XEpMxdXFxlOtqvPiPTXx9lCyT0RujRokEaukJsH6mWZ035ECHsXgnrHQ9ElF97Nw4NQlRi6IYNsfFwBoUNmfqf1CqRfgZ3FlIlJYqEESscIfm2HJYDi73xzX6mJG9/0rW1uXk0tMsTNnzUHeWxNNit3Ax8OVf3WuzaMtquKq+6eJSB5SgyRyO125aJ5ntH2uOfYtdzW631fR/SxsPXSOkQsjOHQmAYB76pRjYp9gKpXQOVoikveKxNXS+vbtS8mSJbn33nuve2zp0qXUrl2bmjVr8vHHH1tQnRQJhgGR38O7zf9qjho/Cs//CsG6rlFmYi+nMOK7XTzw4RYOnUmgbHFP5vyzMR8PaKLmSETyjc0wDMPqIvLb6tWriY+PZ/78+Xz77bdpy1NTU6lXrx6rV6/Gz8+Pxo0bs3XrVkqVyvpaM3Fxcfj7+xMbG4ufn857kEzEHjOvabRvmTkuXQN6vg1V77a2LidnGAZLd51gwpJIzsYnAfBQ8yCGd6mDv7funyYiuZPd7+8icYitXbt2rFmz5rrlv/76K/Xr16dSpUoAdOvWjRUrVvDggw/e5gqlUHI4zKtgr5oAyZfM6P7dL0GrYeDuZXV1Tu3PC5cZs2g3q/edAaBGuWJM7RdC06q6UKaI3B6WH2Jbt24dPXv2JCAgAJvNxqJFi65bZ86cOVSrVg0vLy/CwsJYv359njz38ePH05ojgMqVK3Ps2LE82bcUcaej4NPO5kUfky9B5abw9HpoP1rNUSZS7Q4+Xn+IjjPXsXrfGTxcXXipQy2WDb5bzZGI3FaWzyAlJCTQoEEDHnvsMfr373/d41999RVDhgxhzpw53HXXXXzwwQd07dqVyMhIgoLMm3WGhYWRlJR03bY//fQTAQEBN33uGx1dvNntCJKSktI9R1xcXJavTYqglEQztr/hravR/eLQYRw0eVzR/SzsPhbLyAURRByLBaBZ1VJM6RdCjXLFLK5MRIoiyxukrl270rVr15s+PnPmTJ544gmefPJJAGbNmsWKFSt47733mDp1KgDbt2/P1XNXqlQp3YzRn3/+SfPmzW+47tSpU5kwYUKunkeKiMMbzQs+njtgjmt3g25vgn+lzLcr4i4npzJr1QE+2RCD3WHg5+XGq93qcn+TQFwU3RcRi1h+iC0zycnJbN++nU6dOqVb3qlTJzZt2nTL+2/WrBm7d+/m2LFjXLp0ieXLl9O5c+cbrjty5EhiY2PTfo4ePXrLzy+FxJWLsHgwzOtmNkfFysP9n8E//qPmKAtr95+h01vr+HDdIewOgx6hFVn1chv+0SxIzZGIWMryGaTMnD17FrvdTvny5dMtL1++PCdPnsz2fjp37syOHTtISEigcuXKLFy4kKZNm+Lm5saMGTNo164dDoeDV155hdKlS99wH56ennh66qaX8jeGAZGL4IfhEH/KXBY2EDpMAO8SFhbm/M7GJzFpaSTfhx8HoFIJbyb1qU/7OuWz2FJE5PZw6gbpmoznBRmGcdNzhW5kxYoVN32sV69e9OrVK9e1SREV+ycsGwb7fzDHpWteje7fZW1dTs4wDL7Z/idTlkdx8XIKLjZ47K5qDO1YC1/PAvHXkYgUEU79N1KZMmVwdXW9brbo9OnT180qidwWDjv89jH8PBGS48HFHVoNhVYvg5tmGDMTczaBVxdEsPnQOQDqVfRjWv8QQiuXsLYwEZEbcOoGycPDg7CwMFauXEnfvn3Tlq9cuZLevXtbWJkUSaf2mOcaHdtmjgObm7NG5epaW5eTS0518OG6g8z+JZrkVAde7mZ0//G7q+Hu6tSnQYpIEWZ5gxQfH090dHTaOCYmhvDwcEqVKkVQUBBDhw7lkUceoUmTJrRo0YIPP/yQI0eOMGjQIAurliIlJRHWTYeNs8CRakb3O46HsMfBRV/wmdn+xwVeXRDBvlOXAGhVswxT+oYQWMrH4spERDJneYO0bds22rVrlzYeOnQoAAMGDGDevHk88MADnDt3jokTJ3LixAmCg4NZvnw5VapUsapkKUpi1pvR/fMHzXGdHtBtOvjd/PpaAnGJKUz/cR9fbP0Dw4DSvh6M6VGP3g0DcnT+oIiIVYrEvdjyg+7FVshdPg8rx8LOz81xsQrQ/U2o29PaugqAFXtOMvb73ZyKMy+sel9YZV7tVpeSvh4WVyYionuxieSOYcCeBWZ0P8G8DxhNHocO48HL39LSnN3J2ETGLd7Nij3mJQ+qlvZhSt8QWtYoY3FlIiI5pwZJ5JqLR2HZy3Dg6mUhytQ2T8Ku0sLaupycw2Hw5dY/eP3HfcQnpeLmYmNQmzt4vn0NvNx1exURKZjUIIk47PDrh/DzJEhJAFcPM7Z/90uK7mdh38lLjFywix1HLgLQKKgE0/qFUrtCcWsLExG5RWqQpGg7uRuWDIZjV+/nF9TCnDUqW9vaupxcYoqdd36J5v21B0l1GBTzdGN4l9o81LwKrrpFiIgUAmqQpGhKuQJr34BNs83ovqcfdJwAjQcqup+FTQfPMmrhbmLOJgDQqV55JvYOpoK/l8WViYjkHTVIUvQcWgtLh8D5Q+a4bk/oOh38KlpalrO7kJDM5OVRfLv9TwDK+3kyoVcwXYIrWFyZiEjeU4MkRcfl8/DTGAj/whwXDzCvaVS3h7V1OTnDMPg+/DgTl0ZyPiEZmw0ebl6Ff3WpjZ+Xu9XliYjkCzVIUvgZBuz+zozuXz4L2KDpE3DPOPDSNawyc+TcZUYtimD9gbMA1CpfjKn9QgmrUtLiykRE8pcaJCncLh6BpUMheqU5LlsHes6GoObW1uXkUu0OPtkQw1ur9pOY4sDDzYXB7Wvwf63vwMNN52iJSOGnBkkKJ4cdtr4Pv7wGKZfN6H7rf8FdQ8BNV3TOzK4/LzLiuwgiT8QB0KJ6aSb3DaZ62WIWVyYicvuoQZLC58QuM7p/fKc5Dmp5Nbpfy9q6nFxCUiozftrPvE0xOAzw93ZnVPe63BdWWfdPE5EiRw2SFB7Jl2HtNNj0Dhh28PSHThOh0aOK7mfhl72nGLNoD8cuXgGgd8MAxvSoR5liulCmiBRNapCkcDi4Gpa+BBdizHG93tD1DSiuCHpmTl9KZMKSSJbtOgFA5ZLevNYnmLa1y1lcmYiItdQgScF2+TysGAW//8cc+1WCbm9CnW7W1uXkHA6Dr7cdZcryKOISU3F1sfHE3dUY0qEmPh76a0FERH8TSsFkGBDxDfw4Ai6fA2zQ7CloP0bR/SxEn47n1YUR/BpzHoCQSv5M7RdCcCV/iysTEXEeapCk4Lnwh3k47eDP5rhsXej1bwhsam1dTi4p1c77aw7x7upoku0OvN1deblTLQa2rIqbq87REhH5OzVIUnDYU2Hre7B6ytXovie0+Re0fFHR/Sz8dvg8IxdEEH06HoC2tcvyWp9gKpf0sbgyERHnpAZJCoYTv8PiwXAi3BxXuduM7pepYWlZzi72Sgqv/7iX/2w9AkCZYh6M61mfHqEVFd0XEcmEGiRxbsmXYc0U2DzHjO57+UOn16DRI6Av+JsyDIMfdp9k3OI9nLmUBMCDzQIZ0aUu/j66f5qISFbUIInzOvgLLBkCF/8wx/X7QpfXoXh5S8tydscvXmHs97tZFXUagOplfZnaN4Tm1UtbXJmISMGhBkmcT8I5WPEq7PqfOfarDN1nQO0u1tbl5OwOg882H+bNFftISLbj7mrjmbY1eLbtHXi5u1pdnohIgaIGSZyHYcCur+DHkXDlPGCD5k9D+9HgWdzq6pxa5PE4Ri7Yxe9/xgLQpEpJpvYLoWZ5vW8iIrmhBkmcw/kYM7p/aLU5LlffjO5XDrO2Lid3JdnO2z8f4KP1h7A7DIp7uTGiax0ebBqEi4vO0RIRyS01SGIteypseRdWT4XUK2Z0v+1waDkYXHUycWbWHzjDqIW7OXL+MgDdQiowvmd9yvl5WVyZiEjBpwZJrHN8pxndP7nLHFdtZUb3S99hbV1O7lx8EpOXRbFg5zEAKvp7MbF3MB3r6eR1EZG8ogZJbr/kBPNij1vmgOEArxLQeTI0/Kei+5kwDIMFO47x2rJILlxOwWaDAS2qMqxzbYp56qMsIpKX9Leq3F7Rq8xzjS6aFy4k+F7oMhWK6e7xmTl8NoFRiyLYGH0OgDoVijOtfygNA0tYW5iISCGlBkluj4SzZjot4mtz7B8I3WdCrU7W1uXkUuwOPlp/iLdXHSAp1YGnmwtDOtTiyVbVcNf900RE8o0aJMlfhgG//9e8rtGVC2BzgeaDoN0o8CxmdXVObeeRC4xcEMHek5cAuLtGGSb3DaZKaV+LKxMRKfzUIEn+OX/IvBJ2zFpzXD4Eer0NlRTdz0x8UipvrtjH/M2HMQwo6ePOmB716Nuoku6fJiJym6hBkrxnT4HN78CaaZCaCG5e0HYEtHhe0f0srIw8xdjvd3MiNhGAfo0rMbp7PUr5elhcmYhI0aIGSfLWsR1mdP9UhDmu1gZ6vKXofhZOxSUyfvEefth9EoCgUj5M6RvC3TXLWFyZiEjRpAZJ8kZSvBnd3/qeGd33Lgmdp0CDBxXdz4TDYfCfX4/w+g97uZSUiquLjf9rXZ3B7Wvi7aH7p4mIWEUNkty6Ayth6VCIvRrdD7nfbI6KlbW2Lie3/9QlRi6IYPsfFwBoEFiCaf1CqFvRz+LKREREDZLkXvwZ+HE47P7OHPsHmYfTanawti4nl5hiZ87qaN5be5AUu4Gvhyv/6lybR1pUxVX3TxMRcQpqkCTnDAPCv4QVoyDxohndv/NZaPcqeCiCnpkth87x6oIIDp1NAKBD3XJM7B1MQAlviysTEZG/U4MkOXPuICwdAjHrzHGFUOg1GwIaWVqWs7t4OZmpy/fy1bajAJQr7smEXvXpElxB0X0RESekBkmyx54Cm2bD2jeuRve9od1IuPM5cNWv0c0YhsGSXSeYuGQPZ+OTAfhn8yBe6VIHf29d8kBExFnpm02y9ud2WDIYTu02x9XbQo9ZUKqalVU5vaPnLzPm+92s2XcGgBrlijGtXwhNqpayuDIREcmKGiS5uaRL8Mtk2Po+YIB3KfPGsqEPKLqfiVS7g3mbDjPjp/1cSbHj4erC8+1r8HSb6ni6KbovIlIQqEGSG9u/wozux/1pjkP/AZ0ng68uXJiZ3cdiGbFgF7uPxQHQrFoppvYL4Y6yuu+ciEhBogZJ0rt0yozu71lojktUMaP7Ne6xti4ndzk5lbdW7ueTDTE4DPDzcmNU97rcFxaIi6L7IiIFjhokMRkG7PwcfhoNibFmdL/Fc9B2pKL7WVi97zSjF+7m2MUrAPRsEMDYHvUoW9zT4spERCS31CAJnI02o/uH15vjig2g52wIaGhlVU7vzKUkJi2NZPHvxwGoVMKb1/oE065OOYsrExGRW6UGqShLTYZNb8Pa6WBPAncf82KPzZ9RdD8ThmHwzbY/mbw8itgrKbjY4PG7qvFSx1r4eup9ExEpDPS3eVF19Dczun860hzfcQ/0mAklq1palrM7dCaeVxdGsOXQeQDqB/gxrV8oIZX9La5MRETykhqkoibpEvw8EX79CDDApzR0mQYh9ym6n4nkVAcfrD3Iv1dHk5zqwNvdlaEda/HYXVVxc3WxujwREcljapCKkn0/wLKXIe6YOW7wIHSaDL6lra3LyW3/4zwjvovgwOl4ANrUKstrfYIJLOVjcWUiIpJf1CAVBZdOwg/DIXKROS5Z1bwS9h3tLCzK+cUlpvDGj3v5cusRDANK+3owtmc9ejUI0P3TREQKOTVIhZnDATs/g5/GQlIs2Fyh5QvQZjh4aPYjMz/uPsHY7/dw+lISAPc3qcyr3epSwsfD4spEROR2UINUWJ09AEtehD82muOARmZ0v2KotXU5uROxVxj7/R5WRp4CoFoZXyb3DablHbqCuIhIUaIGqbBJTYaNs2DddLAnm9H99qOh2dOK7mfC7jD4YssfTF+xj/ikVNxcbDzT9g6ea1cDL3fdP01EpKjRN2ZhcvRXWDwYzkSZ4xodoPtMKFnF2rqc3N6TcYz4LoLwoxcBaBxUgqn9Qqldobi1hYmIiGXUIBUGiXHw8wT47RPM6H4Z6Po6BPdXdD8TiSl2Zv98gA/XHSLVYVDc041Xutbhn82CdP80EZEiTg1SQbd3GSwbBpfM213Q8J/Q6TXwKWVtXU5uY/RZRi2M4PC5ywB0qV+B8b3qU8Hfy+LKRETEGahBKqjiTsAPr0DUYnNcshr0fBuqt7G2Lid3ISGZ15ZF8d2OPwGo4OfFhN716Vy/gsWViYiIM1GDVNA4HLBjHqwcb0b3Xdyg5WBo8wq4e1tdndMyDINF4ceYtDSK8wnJ2Gzw6J1VGNa5NsW93K0uT0REnIwapILkzH7z/mlHNpvjgMbQ699QIdjaupzckXOXGbUogvUHzgJQu3xxpvYPoXFQSYsrExERZ1UkbiLVt29fSpYsyb333ptu+dGjR2nbti316tUjNDSUb775xqIKs5CaBGumwft3mc2Ru695/7QnV6k5ykSK3cH7aw/SadZa1h84i4ebC//qXJulg+9WcyQiIpmyGYZhWF1Eflu9ejXx8fHMnz+fb7/9Nm35iRMnOHXqFA0bNuT06dM0btyYffv24evrm+U+4+Li8Pf3JzY2Fj8/v/wr/sgWM7p/dp85rtkJus+AEkH595yFwO9HLzJiQQRRJ+IAaHlHaSb3DaFamaz/bEVEpPDK7vd3kTjE1q5dO9asWXPd8ooVK1KxYkUAypUrR6lSpTh//ny2GqR8lxgLq8bDtk/NsW9ZM7pfv5+i+5mIT0plxk/7mL/pMA4DSvi4M7p7Pfo3rqT7p4mISLZZfoht3bp19OzZk4AA8wagixYtum6dOXPmUK1aNby8vAgLC2P9+vV5Xse2bdtwOBwEBgbm+b5zLGoJvNv8r+ao0cPw3K+6rlEWfo46RaeZa5m70WyO+jaqxM9D23BvWGU1RyIikiOWzyAlJCTQoEEDHnvsMfr373/d41999RVDhgxhzpw53HXXXXzwwQd07dqVyMhIgoLMw0xhYWEkJSVdt+1PP/1EQEBAljWcO3eORx99lI8//vim6yQlJaV7jri4uOy8vJxb9CyEf2n+f6k7oOcsqNY6f56rkDgdl8iEJZEsizgBQGApbyb3CaF1rbIWVyYiIgWV5Q1S165d6dq1600fnzlzJk888QRPPvkkALNmzWLFihW89957TJ06FYDt27fn+vmTkpLo27cvI0eOpGXLljddb+rUqUyYMCHXz5NtAY1g11dw14vQ+l+K7mfC4TD4329HmfpDFJcSU3F1sfFkq2oMuacW3h66f5qIiOSe5Q1SZpKTk9m+fTsjRoxIt7xTp05s2rTplvdvGAYDBw6kffv2PPLII5muO3LkSIYOHZo2jouLy5/DcU2egOptoUzNvN93IRJ9+hIjF0Tw2+ELAIRW9mdqvxDqB/hbXJmIiBQGTt0gnT17FrvdTvny5dMtL1++PCdPnsz2fjp37syOHTtISEigcuXKLFy4kKZNm7Jx40a++uorQkND0859+vzzzwkJCbluH56ennh6et7S68kWFxc1R5lISrUzZ/VB3ltzkGS7Ax8PV4Z1qs2AllVx1f3TREQkjzh1g3RNxhNsDcPI0Um3K1asuOHyu+++G4fDcUu1ye3za8x5Ri7YxcEzCQC0r1OOSX2CqVRChyFFRCRvOXWDVKZMGVxdXa+bLTp9+vR1s0pSeMVeTmHaj1H899ejAJQp5smEXvXpFlJB6TQREckXlsf8M+Ph4UFYWBgrV65Mt3zlypWZnlAthYNhGCzddZx7Zq5Na44ebBbEz0Pb0D20opojERHJN5bPIMXHxxMdHZ02jomJITw8nFKlShEUFMTQoUN55JFHaNKkCS1atODDDz/kyJEjDBo0yMKqJb8du3iFsYt28/Pe0wDcUdaXqf1CaVatlMWViYhIUWB5g7Rt2zbatWuXNr6WFBswYADz5s3jgQce4Ny5c0ycOJETJ04QHBzM8uXLqVKlilUlSz6yOwzmbzrMmz/t43KyHQ9XF55tdwfPtL0DTzdF90VE5PYoEvdiyw+37V5sRcie47GMXBDBrj9jAWhatSRT+4VQo1xxiysTEZHCQvdikwLjSrKdWT/v5+P1MdgdBsW93Hi1W10eaBKIi6L7IiJiATVIYql1+88walEER89fAaB7aEXG9ahHOT8viysTEZGiTA2SWOJcfBKvLYti4c5jAAT4ezGpTzD31NXlG0RExHpqkOS2MgyD73Yc47VlkVy8nIKLDQa2rMbLnWrh66lfRxERcQ76RpLb5vDZBF5dGMGmg+cAqFvRj2n9QmgQWMLawkRERDJQgyT5LsXu4MN1h5j98wGSUh14ubvwUodaPH53NdxdnfpapSIiUkSpQZJ8tePIBUZ+F8G+U5cAaFWzDJP7hBBU2sfiykRERG5ODZLki0uJKby5Yh+fbfkDw4BSvh6M7VGP3g0DdIsQERFxemqQJM/9tOckY7/fw8m4RADuDavMqG51KenrYXFlIiIi2aMGSfLMydhExi/ew497TgJQtbQPk/uGcFeNMhZXJiIikjNqkOSWORwGX/56hDd+2MulpFTcXGw83aY6L7SviZe77p8mIiIFjxokuSX7T11i5IIItv9xAYCGgSWY1j+EOhV0fzoRESm41CBJriSm2Hl3dTTvrz1Iit2gmKcbr3SpzT+bV8FV908TEZECTg2S5Njmg+d4dWEEMWcTAOhYrzwTe9enor+3xZWJiIjkDTVIkm0XLyczZXkUX2/7E4ByxT2Z2DuYLsEVLK5MREQkb6lBkiwZhsHi348zcUkk5xKSsdng4eZV+FeX2vh5uVtdnoiISJ5TgySZOnr+MqMX7Wbt/jMA1CpfjKn9QgirUsriykRERPKPGiS5oVS7g7kbDzNz5X6upNjxcHNhcPsa/F/rO/Bw0/3TRESkcFODJNeJ+DOWEQt2sed4HAB3Vi/FlL4hVC9bzOLKREREbg81SJImISmVt1bu59ONMTgM8Pd2Z1S3utzXpLLunyYiIkWKGiQBYPXe04xetJtjF68A0KtBAGN71qNMMU+LKxMREbn91CAVcWcuJTFxaSRLfj8OQOWS3rzWJ5i2tctZXJmIiIh11CAVUYZh8PW2o0xeFkVcYiouNnji7mq81LEWPh76tRARkaJN34RF0MEz8by6IIKtMecBCK7kx7R+oQRX8re4MhEREeegBqkISU518P7ag7zzSzTJdgfe7q683KkWA1tWxc1V0X0REZFr1CAVEdsOn2fkgggOnI4HoG3tskzqHUxgKR+LKxMREXE+apAKudgrKbzx416+3HoEgDLFPBjXsz49Qisqui8iInITapAKKcMw+HH3ScYt3sPpS0kA/KNpICO61qGEj4fF1YmIiDg3NUiF0PGLVxj7/R5WRZ0CoHoZX6b0C+HO6qUtrkxERKRgUINUiNgdBp9vPsz0FftISLbj7mrjmTZ38Gy7Gni5u1pdnoiISIGhBqmQiDoRx4gFEfx+9CIATaqUZGq/EGqWL25tYSIiIgWQGqQCLjHFzts/H+CjdYdIdRgU93RjeNc6PNQsCBcXnYQtIiKSG2qQCrANB84yalEEf5y7DEDX4AqM71Wf8n5eFlcmIiJSsKlBKoDOJyTz2rJIFuw4BkAFPy8m9QmmY73yFlcmIiJSOKhBKkAMw2DhzmNMWhrJhcsp2GwwoEVVXu5Ui+Je7laXJyIiUmioQSog/jiXwKiFu9kQfRaAOhWKM7VfCI2CSlpcmYiISOGjBsnJpdgdfLw+hlmr9pOU6sDTzYUXO9TkqVbVcdf900RERPKFGiQnFn70IiO+28Xek5cAuKtGaSb3CaFqGV+LKxMRESnc1CA5ofikVN5csY/5mw9jGFDSx53R3evRr3El3T9NRETkNlCD5GRWRZ5izPe7ORGbCEC/RpUY1b0upYt5WlyZiIhI0aEGycl8//txTsQmElTKh8l9g2lVs6zVJYmIiBQ5apCczNge9aha2odn29bA20P3TxMREbGCGiQnU7a4Jy93qm11GSIiIkWacuIiIiIiGahBEhEREclADZKIiIhIBmqQRERERDJQgyQiIiKSgRokERERkQzUIImIiIhkoAZJREREJAM1SCIiIiIZqEESERERyUANkoiIiEgGapBEREREMlCDJCIiIpKBm9UFFFSGYQAQFxdncSUiIiKSXde+t699j9+MGqRcunTpEgCBgYEWVyIiIiI5denSJfz9/W/6uM3IqoWSG3I4HBw/fpz27duzbdu2bG3TtGlTfvvtt0zXiYuLIzAwkKNHj+Ln55cXpRZ42XnfrHS768uv58ur/d7KfnKzbU630ecwd5z5c2hFbc78Obzdn8Gcbmf1Z9AwDC5dukRAQAAuLjc/00gzSLnk4uJC5cqVcXNzy/Yfnqura7bX9fPz01/MV+XkfbPC7a4vv54vr/Z7K/vJzbY53Uafw9xx5s+hFbU58+fwdn8Gc7qdM3wGM5s5ukYnad+i5557Ll/Wlb84+/t2u+vLr+fLq/3eyn5ys21Ot3H23ydn5czvmxW1OfPn8HZ/BnO6nTP/Lv2dDrE5mbi4OPz9/YmNjXXaf62JFHb6HIpYyxk+g5pBcjKenp6MGzcOT09Pq0sRKbL0ORSxljN8BjWDJCIiIpKBZpBEREREMlCDJCIiIpKBGiQRERGRDNQgiYiIiGSgBklEREQkAzVIBdilS5do2rQpDRs2JCQkhI8++sjqkkSKlKNHj9K2bVvq1atHaGgo33zzjdUliRRJffv2pWTJktx77715tk/F/Aswu91OUlISPj4+XL58meDgYH777TdKly5tdWkiRcKJEyc4deoUDRs25PTp0zRu3Jh9+/bh6+trdWkiRcrq1auJj49n/vz5fPvtt3myT80gFWCurq74+PgAkJiYiN1uR/2uyO1TsWJFGjZsCEC5cuUoVaoU58+ft7YokSKoXbt2FC9ePE/3qQYpH61bt46ePXsSEBCAzWZj0aJF160zZ84cqlWrhpeXF2FhYaxfvz5Hz3Hx4kUaNGhA5cqVeeWVVyhTpkweVS9S8N2Oz+A127Ztw+FwEBgYeItVixQut/NzmJfUIOWjhIQEGjRowDvvvHPDx7/66iuGDBnCqFGj2LlzJ61ataJr164cOXIkbZ2wsDCCg4Ov+zl+/DgAJUqU4PfffycmJob//Oc/nDp16ra8NpGC4HZ8BgHOnTvHo48+yocffpjvr0mkoLldn8M8Z8htARgLFy5Mt6xZs2bGoEGD0i2rU6eOMWLEiFw9x6BBg4yvv/46tyWKFGr59RlMTEw0WrVqZXz22Wd5UaZIoZaf34WrV682+vfvf6slptEMkkWSk5PZvn07nTp1Sre8U6dObNq0KVv7OHXqFHFxcYB55+N169ZRu3btPK9VpDDKi8+gYRgMHDiQ9u3b88gjj+RHmSKFWl58DvOLm6XPXoSdPXsWu91O+fLl0y0vX748J0+ezNY+/vzzT5544gkMw8AwDJ5//nlCQ0Pzo1yRQicvPoMbN27kq6++IjQ0NO28is8//5yQkJC8LlekUMqLzyFA586d2bFjBwkJCVSuXJmFCxfStGnTW6pNDZLFbDZburFhGNctu5mwsDDCw8PzoSqRouNWPoN33303DocjP8oSKVJu5XMIsGLFirwuSSdpW6VMmTK4urpe1yGfPn36uk5aRPKePoMi1nPmz6EaJIt4eHgQFhbGypUr0y1fuXIlLVu2tKgqkaJDn0ER6znz51CH2PJRfHw80dHRaeOYmBjCw8MpVaoUQUFBDB06lEceeYQmTZrQokULPvzwQ44cOcKgQYMsrFqk8NBnUMR6BfZzmGd5OLnO6tWrDeC6nwEDBqSt8+677xpVqlQxPDw8jMaNGxtr1661rmCRQkafQRHrFdTPoe7FJiIiIpKBzkESERERyUANkoiIiEgGapBEREREMlCDJCIiIpKBGiQRERGRDNQgiYiIiGSgBklEREQkAzVIIiIiIhmoQRIRERHJQA2SiBQ6bdu2ZciQIfn6HOPHj6dhw4b5+hwiYh01SCIieWDgwIH06dPH6jJEJI+oQRIRERHJQA2SiBRqFy5c4NFHH6VkyZL4+PjQtWtXDhw4kG6djz76iMDAQHx8fOjbty8zZ86kRIkS6daZNm0a5cuXp3jx4jzxxBMkJiamPTZ+/Hjmz5/P999/j81mw2azsWbNmtvw6kQkv6hBEpFCbeDAgWzbto3FixezefNmDMOgW7dupKSkALBx40YGDRrEiy++SHh4OB07dmTy5Mnp9vH1118zbtw4Jk+ezLZt26hYsSJz5sxJe3zYsGHcf//9dOnShRMnTnDixAlatmx5W1+niOQtm2EYhtVFiIjkpbZt29KwYUOee+45atWqxcaNG9MalnPnzhEYGMj8+fO57777+Mc//kF8fDxLly5N2/7hhx9m6dKlXLx4EYCWLVvSoEED3nvvvbR17rzzThITEwkPDwfMRuzixYssWrTodr1MEclHmkESkUIrKioKNzc3mjdvnrasdOnS1K5dm6ioKAD27dtHs2bN0m2XcRwVFUWLFi3SLcs4FpHCRQ2SiBRaN5sgNwwDm8123f9ntZ2IFB1qkESk0KpXrx6pqals3bo1bdm5c+fYv38/devWBaBOnTr8+uuv6bbbtm1bunHdunXZsmVLumUZxx4eHtjt9rwsX0QspAZJRAqtmjVr0rt3b5566ik2bNjA77//zsMPP0ylSpXo3bs3AC+88ALLly9n5syZHDhwgA8++IAffvgh3azSiy++yKeffsqnn37K/v37GTduHHv27En3XFWrVmXXrl3s27ePs2fPpp0ELiIFkxokESnU5s6dS1hYGD169KBFixYYhsHy5ctxd3cH4K677uL9999n5syZNGjQgB9//JGXXnoJLy+vtH088MADjB07luHDhxMWFsYff/zBM888k+55nnrqKWrXrk2TJk0oW7YsGzduvK2vU0TyllJsIiIZPPXUU+zdu5f169dbXYqIWMTN6gJERKz25ptv0rFjR3x9ffnhhx+YP39+uusciUjRoxkkESny7r//ftasWcOlS5eoXr06L7zwAoMGDbK6LBGxkBokERERkQx0kraIiIhIBmqQRERERDJQgyQiIiKSgRokERERkQzUIImIiIhkoAZJREREJAM1SCIiIiIZqEESERERyeD/Ab8z2kdViLawAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#plotting error vs array dt_array\n",
    "\n",
    "def error_plot(dt_array):\n",
    "\n",
    "    vmy=np.vectorize(numerical_solution_error)\n",
    "\n",
    "    vlib=np.vectorize(library_error)\n",
    "\n",
    "    myerror=vmy(dt_array) #error of my rk-4 code\n",
    "\n",
    "    liberror=vlib(dt_array) #error from library function\n",
    "\n",
    "    #plotting my code error\n",
    "\n",
    "    plt.loglog(dt_array,myerror,label='actual error')\n",
    "\n",
    "    #plotting the expected rk-4 slope\n",
    "\n",
    "    plt.loglog(dt_array,dt_array**4,label='expected slope')\n",
    "    \n",
    "    #plotting library code error\n",
    "\n",
    "    plt.loglog(dt_array,liberror,label='library error')\n",
    "\n",
    "    plt.legend()\n",
    "    plt.xlabel('logdt')\n",
    "    plt.ylabel('logerror')\n",
    "    plt.title('error plot')\n",
    "\n",
    "dt_array=np.array([0.1,0.01,0.001])\n",
    "\n",
    "error_plot(dt_array)\n",
    "\n",
    "\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "#the odeint doesnt feel like rk4 method \n",
    "\n",
    "#couldnt find a rk4 method which depends on stepsize .checked rk45 but the step size is chosemn by solver"
   ]
  }
 ],
 "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}