m-thesis-introduction/fourier/03_phase_by_correlation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Phase by Correlation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Continuous Cross Correlation\n",
    "For signals $u(t)$ and $v(t)$ \n",
    "and correlation parameter $\\tau$:\n",
    "\n",
    "$$\n",
    " Corr(\\tau; u, v)\n",
    "  =\n",
    "  \\int_{-\\infty}^\\infty \\mathrm{d}t \\; u(t) \\, v^*(t - \\tau)\n",
    "$$\n",
    "\n",
    "\n",
    "##### Discrete Cross Correlation\n",
    "For signals $u[n]$ and $v[n]$ \n",
    "and correlation parameter $k$:\n",
    "\n",
    "$$\n",
    " Corr(k; u, v)\n",
    "  = \\sum_n u[n] \\, v^*[n-k]\n",
    "$$\n",
    "\n",
    "---\n",
    "\n",
    "Writing\n",
    "\n",
    "$$\n",
    " u(t, \\phi_t) = A_1 \\exp(2i \\pi f t + \\phi_t)\n",
    "$$\n",
    "\n",
    "\n",
    "Correlation $u \\triangleq u(t_1, \\phi_{t1})$ and $ v \\triangleq u(t_2, \\phi_{t2}) $:\n",
    "\n",
    "$$\n",
    " Corr(\\vec{\\tau}; u_1, u_2)\n",
    "  =\n",
    "  \\int_{-\\infty}^\\infty \\mathrm{d}\\vec{\\zeta} \\; u_1(\\vec{\\zeta}) \\, u_2^*(\\vec{\\zeta} - \\vec{\\tau})\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import signal\n",
    "rng = np.random.default_rng()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "class SineWave:\n",
    "    \"\"\"\n",
    "        A periodic sine wave, totally described by its frequencies, phases and amplitudes.\n",
    "    \"\"\"\n",
    "    def __init__(self, frequencies, phase = 0, amplitude = 1):\n",
    "        self.freqs = frequencies\n",
    "        self.phase = phase\n",
    "        self.amplitude = amplitude\n",
    "    \n",
    "    def __call__(self, time, complx=True):\n",
    "        if complx:\n",
    "            return self.amplitude * np.exp( 2j * self.freqs*time + 2j * self.phase)\n",
    "        else:\n",
    "            return self.amplitude * np.cos( self.freqs*time + self.phase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Initial signal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "mysignal = SineWave( 5e1, np.pi/2, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Frequencies: 50.0\n",
      "Phase: 0.5π\n",
      "Amplitude: 1\n"
     ]
    }
   ],
   "source": [
    "# show the signals\n",
    "print(\"Frequencies:\", signal.freqs)\n",
    "print(\"Phase:\", \"{}π\".format(signal.phase/np.pi))\n",
    "print(\"Amplitude:\", signal.amplitude)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x14591292bb10>]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOy9d3hUx9m/f8+qFyQQaiABoghUEFVUYxtTTTHgXmM7ceIkb5z+S+LkTU+cOMk3id/0OHYc3AvGBoMwohiMbZooQh0JUSTUJYSEunbn98esHBlLqO3u2ZHOfV177e7ZUz6I3fPMPPMUIaXExMTExGToYjFagImJiYmJsZiGwMTExGSIYxoCExMTkyGOaQhMTExMhjimITAxMTEZ4ngaLaA/hIaGypiYGKNlmJiYmGjFsWPHqqSUYVdv19IQxMTEkJaWZrQMExMTE60QQpzvarvpGjIxMTEZ4piGwMTExGSIYxoCExMTkyGOaQhMTExMhjimITAxMTEZ4jjEEAgh/i2EqBBCZHbzuRBC/EkIUSCEOCWEmNXps4eEEPn2x0OO0GNiYmJi0nscNSP4D3DzNT5fBcTaH48CfwcQQoQAPwHmAXOBnwghRjhIk4mJiYlJL3BIHoGU8n0hRMw1dlkPPC9VzetDQojhQohRwGJgl5SyBkAIsQtlUF5xhK6r+c2R35Bbk+uMU5t0g01KrjS309hmxWqTCMDXy4Nhvp54eWjgmbS1Q0s9tDWCzQYWD/DyA58g9drEpBsk0NDSTmOrlTarDQAfTwsBPp74efXvuxMXEsf35n7PgSoVrkooiwKKOr0vtm/rbvunEEI8ippNMHbsWOeoNHEYTa1WSmqbqGlsxWr7dM8LISDQx4tRwb6EBHgboLAHmmvhcgk0X4KuenYIAX4hEBQFvkGu12fitrRZbZRebqbqSgut7bYu9/Hz9iBimC/hQT5YhHCxwk/jKkPQ1b9UXmP7pzdK+TTwNEBycnK/uuk4w5KafJLLTW38OiWHV48WEeDtwboZUaxIiCBhdBAj/L1pbrdypuIK75+uYvOJYtJPN5IUFcyTtyeRODrYaPlw6Txs+wac2QuBkTD9XohdDuEJ4DMMmmqhIhvyU+Hky1CYA/G3wM1PQnC00epNDMRqkzz7QSF/3JVPS7uVZfERrJk2ihljhhMR5IuUcKGmkSNnq9lysoS0E5eoHubDz9cncvPUUYZqF47qUGZ3DW2TUk7t4rN/AvuklK/Y3+eh3EKLgcVSyi92tV93JCcnS7PEhPuRXlTLl148RnldM48sGs//LJ7EiGuM9tutNt45VcIT23OpbWzlB6vj+ex1MQijRkiZb8LWr6vXN30fkh8BL9/u929thMN/h/2/A08fuO1pmLzSNVpN3IqK+mb+58XjpJ2/xLL4CH6wOo4JYYHXPObI2Rp+ujWL7NI67pgdzS83TMW3ny6j3iKEOCalTP7UdhcZgjXAY8Bq1MLwn6SUc+2LxceAjiii48DsjjWD7jANgfux7VQJ33o9nbBAH/52/yymjxne62NrG1v57qZTpGaXc1dyNL+6NQlPV64fSAn7fwv7fgVj5qsb+ohxvT+++gy88RCUZaiZwfwvO0+riduRW1bHZ587Sm1jG7+6bSobZkT1ejDTZrXxpz35/HlvAdOjg3nus3Od6irtzhA4Knz0FeAgMEUIUSyEeEQI8SUhxJfsu6QAhUAB8C/gfwDsN/xfAEftj5/3ZARM3I9Nx4r52isnmB4dzDtfXdQnIwAw3N+bfzwwm68umcTracV8/dWTHy+uOR0pIfWHyghMvw8e2to3IwAwciI8sku5iN59HPY96RytJm5H5sXL3PP0IaSETV9ewK0zo/s0o/XysPDtFVP4xwOzySmr556nD1JR3+xExV3jsBmBKzFnBO5DSkYpX3n5ONdNDOVfDybj5z2wqe2/3i/kiZQcbpsZxe/vmu58N9GeX8CB/wdzvgCrf6cWgfuLzQpbvwonX4IVT8DCxxyn08TtKKio5/a/HyTQx5NXvjCfsSP9B3S+j85U8fmNaYwbGcBrX5xPkK+Xg5T+F6fOCEyGJkfO1vCN104ya+wInnlo4EYA4As3TODbyyez+cRFfrszzwEqr8Hx55URmPXgwI0AqHDSdX+GhA2Q+r9qzcFkUFJR38xD/z6Kl4dwiBEAWDgxlH9+Zjb55fU8+nxatxFHzsA0BCb9ovRyE19+8RjRw/145sFkhy5yPbZkEvfOHcvf951ha3qJw877Cc4fhG3fhIlLYc0fB24EOrB4qDWGsQvg7a9AyUnHnNfEbWiz2vifF49T09DKcw/PdYgR6OD62DB+d+c0DhXW8MT2bIedtydMQ2DSZ1rbbXz5xeM0t1l5+sHZ14wM6g9CCH62LpHkcSP43qZT5JXVO/T8XKmETZ+F4ePgzufAw8FR1J4+cNfz4B8Cr38Gmi879vwmhvKrlBzSzl/iN3dMIyna8SHPt86M5vOLxrPx4HneOlHs8PN3hWkITPrMU7tPc7Kolt/eMZ1J4cOccg1vTwt/u38WAT4efP3VEzS3WR1zYinhrS9C0yV1s/Z1Uu5CYDjcuREuX4R3vtF1UpqJduzNLee5D8/x8MIY1k0f7bTrPL4qjrnjQ/jhW5mcr25w2nU6MA2BSZ84craGv+8/w93JY1gzzblJMOFBvvzujunkltXzO0etF6Q9C2f2wIpfQuSnIp0dy5g5cNMPIGsznHrdudcycTrVV1r47qYM4iKH8f3VcU69lqeHhT/ePQOLRfCN107S7uQoOtMQmPSa5jYr39mUztgQf358S4JLrnlTXDgPzB/Lvz88y7HzlwZ2spqzkPojtS4w5/OOEdgTi74JY+bBju9CfblrrmniFH6yNYu6pjaeumcGPp7OrzMVNdyPX26YyokLtTz34TmnXss0BCa95qnd+ZyvbuTXtyUR4OOq6iTw+Kp4RgX58vibp/ofSSElbP8WCHtkj6uyly0esO4v0NYEO77jmmuaOJw9OeVsO1XKV5dMIi7SdbWl1k0fzbL4cH6/K48L1Y1Ou45pCEx6RV5ZPf86UMhdydEsnBjq0msH+njyy1unkl9xhX8dKOzfSTLfVPWDlv4Igrusa+g8wibDjd+B7C1QuM+11zYZME2tVn70diaTIwL54o0TXXptIQS/2DAVT4uFH23JxFl5X6YhMOkVT6TkEOjjyfdXxRty/SVxEaxIiOBv7xX0PfOytVFlD4+e6TqX0NUs+KqKUnr3+2BtN0aDSb/414FCSi4388sNSXh7uv6WOSrYj28un8z+05XszqlwyjVMQ2DSI+/lVfD+6Uq+tjTW4aGifeH7q+Nptdr4Q+rpvh148K9QX6rqABnVQ8DLVy1QV2TD8f8Yo8Gkz5TXNfP3fWdYNTWSueNDDNPx4IJxxIYH8ott2Y6LoOuEaQhMrkm71cYT23OIGenPZ+b3sQaPgxkfGsCDC2J4La2I7JK63h1UXw4f/BHi1sLY+c4V2BPxt0DM9bD3CWg0S2rpwO9T82i32Xh8lXOjhHrCy8PCT9cl0thq5UzlFYef3zQEJtfklaNFFFRc4fur4w2ZFl/N15bEEuznxS+3Z/fOX7rv12BtgWU/c764nhACbv61anqz/zdGqzHpgeySOt44VszDC2MYNzLAaDlcNymUA9+9ySl9O4z/ZZu4LXXNbfxx12nmjQ9hRUKE0XIACPb34htLY/noTDV7c3vwl1bkqnpCyZ+D0EmuEdgTkUmqttHRZ6D2gtFqTLpBSskTKdkE+3nx2E2xRsv5GEfU8+oK0xCYdMt/PjxHTUMrP1yTYFyzmC64f/44Ykb684ddp689K9jzM/AOgBvdrDPdDd8BYYEDvzdaiUk3HCys5sOCajUD9Xd8FVB3wzQEJl1S39zGsx+cZVl8uFPqqQwELw8LX7lpElkldd3PCkrTIS8FFn4VAlwb7tojwdFqVnDiRXNW4Kb8eU8B4cN8uG/e0OiPbhoCky554dB5Lje18dUl7jMt7syGmVGMCfHjT3vyu54VHPg9+ATB3EddL643LPqmOStwU46eq+FgYTWP3jDB6a0j3QVHdSi7WQiRJ4QoEEI83sXnfxRCnLQ/Tgshajt9Zu302VZH6DEZGI2t7Txz4CyLp4T1uduYq/DysPCVxZNIL77M/tOVn/ywMg+yt8LcL4Cfe+o3ZwXuy5/25BMa6M3984yNknMlAzYEQggP4K/AKiABuFcI8YlCNFLKb0opZ0gpZwB/BjZ3+rip4zMp5bqB6jEZOC8dukBNQ6vbzgY6uG1WNFHDu5gVHPgDePnB/P8xTlxv+HhW8AejlZjYOVlUy4H8Kj5//QSnLcy6I46YEcwFCqSUhVLKVuBVYP019r8XeMUB1zVxAs1tVv75fiGLJoUye9wIo+VcE29PC19aPJHjF2r5sKBabaw5CxlvwOzPut/awNUER8PMz6hZQV2p0WpMgD/vyWe4vxcPGJwz42ocYQiigKJO74vt2z6FEGIcMB7Y22mzrxAiTQhxSAixobuLCCEete+XVllZ2d1uJgNk07Fiqq608NgSNwm37IG7kqOJCPLhn++fURs+fEplDy/8qrHCesvCx8DWDkf/ZbSSIU9uWR17cit45LrxBLqwqKI74AhD0FVcYXcxffcAm6SUnXOkx9qbKd8HPCWE6LKqk5TyaSllspQyOSwsbGCKTbrEZpP8+8OzTIsOZp6B6fR9wcfTgwcXxHAgv4ozZ8/CyZdhxv0Q5NxeCQ4jZALErYG0f0Or8xuQmHTPcx+cw9fLwmcWDK3ZADjGEBQDYzq9jwa6azR7D1e5haSUJfbnQmAfMNMBmkz6wf78SgorG3hk0Xi3yhvoifvmjsXH00Lhu38Ga6v7rw1czYLHVMe0dNNjahTVV1p46+RFbp8VzXB/4+ppGYUjDMFRIFYIMV4I4Y262X8q+kcIMQUYARzstG2EEMLH/joUuA5wXcdmk0/w7w/OEhHkw6qpmoym7YwI8ObOGRFML9tMa8xNquyzToydD6NnwcG/gc25nahMuublwxdobbfx2etijJZiCAM2BFLKduAxYCeQA7wupcwSQvxcCNE5Cuhe4FX5yaDveCBNCJEOvAc8KaU0DYEBnC6v50B+FQ8uiHGLmkJ95bGITMLFJXYGXitOwU0RAhZ8BWrOQP5Oo9UMOVrbbTx/6Dw3Tg5zWg9ud8chKyJSyhQg5aptP77q/U+7OO4jIMkRGkwGxnMfnsPH08K9c/XMpIzM/Q+lntH8MjeKle02/YxZwnrY9RNVMnvKKqPVDCm2Z5RQWd/C5+4cb7QUw9Ds12LiDC41tLL5eDG3zYoixMB+A/2m6ChcPMaV6Z+j/Eob2051t0Tlxnh4wbxH4dwBKM8yWs2QQUrJsx+cZVJ4IDfEunm4sRMxDYEJr6UV0dJu4+GFmo6IDv8DfIKYtPwLTAwL4IVD541W1D9mPAAe3nDsP0YrGTIcv1BL5sU6Hl4Yo1WAhKMxDcEQx2aTvHrkAnNjQpgSqaF/9EoFZL8NMx9A+AZx37xxnLhQS05pLxvXuBMBIyFhA6S/aoaSuohXjlwgwNuDW2e6uI+1m2EagiHOocJqzlU3cu+8MT3v7I6cfFklZM1+GIDbZ0Xh7WnhlSOa1u9J/hy01EHm5p73NRkQl5uUG3H9zCgChlgC2dWYhmCI8/KRCwT7eWkXMgqAlKrxzNgFEDYFgOH+3qxJGsVbxy/S1Or43q5OZ+x8CItTCWYmTmXLyYs0t9m4T9MACUdiGoIhTNWVFnZmlXH7rGg9y+2e+0CFXM566BOb7507lvqWdt7RcdFYCDUrKDkOJSeNVjNokVLy8uELJEUFMzXKvfptGIFpCIYwbx4rps0quU9Xt9Cx/4BvMCR+skTVnJgRTAoP1Nc9NO1u8PSDY88ZrWTQcqKoltyy+iHTeKYnTEMwRJFS8op9kVjLJJrGGsjZqm6aXn6f+EgIwb1zx+q7aOw3HKbeDqfegJZ6o9UMSl45rBaJb5k+2mgpboFpCIYoB3VfJE5/VdUVusot1IH2i8azH4K2BsjeYrSSQUddcxvvnCph3YyoIVdltDtMQzBEeSOtmGG+nhovEm+EqNkQObXLXYb7e3NzYiRbTpbQ0q7honH0HBg5SUVFmTiUbemlNLfZuGeOpoMgJ2AagiHIlZZ23s0sY+200XouEpcch8pc1erxGtw+O5rLTW3szemmwb07IwRMvxfOf6ia7Zg4jDePFxMbHsi0aHORuAPTEAxBUjJKaWqzcsdsTZNo0l8FDx9IvPWauy2aFEpEkA9vHi92kTAHM/0eQKh/r4lDOFvVwLHzl7h9dvSQziS+GtMQDEHePFbM+NAAZo1171aUXdLeChmbVDMX32uP6Dwsgg0zo3gvr5LK+hYXCXQgwdEwYTGkv2yWp3YQm48XYxEM+UziqzENwRCjqKaRw2druG1mlJ4jooJd0FSj3Ca94I5Z0Vhtki0nLzpZmJOYcT/UXlAuIpMBYbNJNh+/yKLYMCKCfI2W03fqSuDD/1MRcw7GNARDjM3H1Q3x1lmajojSX4GAMJi4pFe7x0YMY1p0MG8e19QQxK0BnyBz0dgBHD5bw8XaJm7X9bt/6nXY9WPVzc7BmIZgCCGlZPOJYhZMGEn0CH+j5fSdxhrIexeS7gKP3of93T4rmpzSOrJLNMwp8PZXCXPZW6DlitFqtObN48UE+niyIiHSaCn949TrED0XRnbZ1n1AOMQQCCFuFkLkCSEKhBCPd/H5w0KISiHESfvj850+e0gIkW9/dB0UbuIQ0s5f4nx1I7fPjjZaSv/I2gy2Nvsiau+5ZfpovDyExovG96mcgryUnvc16ZLG1nZ2ZJSyJmkUft4aRsqVZUBFFky7yymnH7AhEEJ4AH8FVgEJwL1CiIQudn1NSjnD/njGfmwI8BNgHjAX+IkQQsMVTD1468RF/Lw8uHmqpiOi9FchPBEi+9bULiTAm5umhLM1vQSrTfZ8gLsxZh4Ej1EjQpN+kZpVTkOrldu0dQu9BhZPlXHuBBwxI5gLFEgpC6WUrcCrQG8bx64Edkkpa6SUl4BdwM0O0GRyFW1WGzsySlmWEKFnNmX1GSg+CtPvVjH2fWT9jCgq61s4VFjtBHFOxmJRN4Aze6Ghymg1WrI1vYTRwb7MiQkxWkrfsVlVpFzsCvB3jn5HGIIooKjT+2L7tqu5XQhxSgixSQjRkdLX22MRQjwqhEgTQqRVVlY6QPbQ4oP8Ki41trFO19oqHfX5+zkiWhofToC3B1tPaliRFCDpTpBWyHrLaCXacamhlfdPV3LL9NFYLBpGyp3/COpLIekOp13CEYagq7/s1fPvd4AYKeU0YDewsQ/Hqo1SPi2lTJZSJoeFhfVb7FBly8mLBPt5ceNkDf92UkLmJhi7UMXW9wNfLw9WJkaSklmqZ8mJiEQIi1cjQ5M+kZJZSrtNsm6GroOgN8ErACY7z1niCENQDHQu2hENfGLYJaWsllJ2ZPT8C5jd22NNBk5Tq5XU7HJWTY3E21PDQLHyLFVSImlg/tFbZoymvrmd/XkaziiFgGl3QtEhuKRpT2aD2HqyhIlhASSMCjJaSt+xtqmIsSmrwDvAaZdxxF3hKBArhBgvhPAG7gG2dt5BCNG5stk6IMf+eiewQggxwr5IvMK+zcSB7Mktp7HVqvGIaBMID9XPdwAsmhTKCH8vtqZrOtbocItlvmmsDo0ovdzEkXM1rJuuaQLlmfdUAqUT3ULgAEMgpWwHHkPdwHOA16WUWUKInwsh1tl3+5oQIksIkQ58DXjYfmwN8AuUMTkK/Ny+zcSBbDlZQvgwH+aNH2m0lL4jpbrxTVgMAaEDOpWXh4U100axO6echpZ2h8hzKSNiVASR6R7qNdvSS5ESjQdBb4LvcJi41KmXcYifQEqZIqWcLKWcKKV8wr7tx1LKrfbX35dSJkopp0spb5JS5nY69t9Sykn2h9mSycFcbmpjf14la6eNxkPHhbLiNFViwUEjonXTo2hus7Eru9wh53M5SXeqePLyLKOVaMHW9BKSooIZH+o8t4rTaGuC3G2QsA48vZ16KQ0dxiZ9YWdmGa1WG+u1HRFtUpVG49Y45HTJ40YwOthX39pDCetBWCDrbaOVuD2FlVfIuHhZ30i50zuh9QpMda5bCExDMOjZllHK2BB/PWuv2+zhkrHLe6w02lssFsHa6aP5oKCKy01tDjmnSwkMh3HXQfbbym1m0i3bT5UCsHa6hs2XQLmFAiMgZpHTL2UagkFMbWMrHxVUsTpplJ4LZec/givlDs+mXDU1kjarZLeu7qGE9VB1Gipyet53CJOSWcbscSMYFezX887uRks95KeqAAmL80timIZgEJOaXU67TbI6SdOSEtlbwNNPZVQ6kBljhjM62JeUjFKHntdlxK8DhNnP+BqcrWogp7SOVbqWU8lPhfbmHpsvOQrTEAxiUjJKiR7hR1KUjm4hG+Rshdhl4BPo0FMLIViVNIoD+VXUNWvoHhoW8V/3kEmXdBj5VUmauoWyt0BgpIoScwGmIRikXG5s40Od3UJFh5VbaIC5A92xOmkUrVYbe3I0dQ8lblBJdhW5Pe87BNmRWcqMMcOJGq6hW6i1AfJ3Qfwtqs6UCzANwSBld045bVap79Q4e4uKFnKwW6iDmWOGExnkS0pGmVPO73Tib0G5h8xZwdVcqG4k82Ida3SdDRTshrZGtRbkIkxDMEhJyShldLAvM8YMN1pK3+lwC01aCr7OKQtgsQhWJUWy/3Ql9Vq6hyJh7AIzjLQLUjKVW0jbcuvZW8A/FMYtdNklTUMwCKlrbuNAfhWrdHULXTwGdRedPiJanTSK1nYbe3MrnHodp5G4ASpzoDLPaCVuRUpGKdOigxkTomEXvrYmlT8Qv9Yl0UIdmIZgELInp5xWq43Vuk6Ns98Gi5cqtOVEZo8dQfgwH/2jh8xZwccU1TRyqviyvt/9M3tVEpkL3UJgGoJBSUpGGZFBvszU0S0kJWRvVc3pHZRE1h0Wi2DV1Ej25VXqWXsoaBSMnW+uE3Ti3Uy15rN6qqaGIHsL+I2AmOtdelnTEAwy6pvb2H+6klVJkXo24Sg5AZcvuGxEtDppFC06u4cSNkBFNlSeNlqJW7A9o5SpUUGMHamhW6i9BfJ2qHIqHl4uvbRpCAYZe3MraG3X2S20RfVmdbJbqIPkmBBCA33YkampeyjBXuDXTC7jYm0TJ4tqWaXrbKBwH7TUOS1k+lqYhmCQkZJRSvgwH2aPHWG0lL4jpbqhjb/Rab1Zr8bD7h7am1tBY6uO7qHRMGa+aQiAHfa1Hq0HQT7B6vvvYkxDMIhoaGlnX14lq6Zq6hYqy4BLZ12+ULYqKZLmNhvv5WrYuQxUTkF5Blw6Z7QSQ9mRWUb8qCA9S063t6qS03GrnV5yuiscYgiEEDcLIfKEEAVCiMe7+PxbQohse/P6PUKIcZ0+swohTtofW68+1qT3vH+6kpZ2GzfrOjXO3qI6kcWtdell540fSUiAN6nZmiaXdZToztlmrA4Dqahr5tj5S/omUJ57H5ovu3wQ1MGADYEQwgP4K7AKSADuFUIkXLXbCSDZ3rx+E/DbTp81SSln2B/rMOk3qdnlDPf3Yk6Mhm4hgJx3IOY6CHBtJzUPi2BpXPjH6yvaETIeIqZC7najlRjGLnupkJWJmhqC7K3gHQgTbjLk8o6YEcwFCqSUhVLKVuBV4BNmTUr5npSy0f72EKpJvYkDabOqyJelcRF4emjo8avKh6o8iLvFkMuvTIykvrmdQ4XVhlx/wMSthQsH4Yqm7q0BsjOrnHEj/Zkc4dgChS7BZlPRQrHLwcvXEAmOuGNEAUWd3hfbt3XHI8COTu99hRBpQohDQohul8uFEI/a90urrByaX/ZrcfRsDZeb2lieEGG0lP7RMZqNW23I5RfFhuLv7aGveyh+LSAhL8VoJS6nrrmNg2eqWJkYqWkmfRo0VLjcJdoZRxiCrv7yXbZOEkI8ACQDv+u0eayUMhm4D3hKCDGxq2OllE9LKZOllMlhYWED1TzoSM0ux8fTwg2TB9bg3TByt8Oo6RBszGTR18uDGyeHkZpVjs2mYeeviKkwfOyQdA/ty6ukzSpZoe0gaJvKpI9dbpgERxiCYmBMp/fRQMnVOwkhlgH/C6yTUrZ0bJdSltifC4F9wEwHaBpSSCnZlV3O9bFh+Ht7Gi2n79SXQfFRQ0dEACsSI6iobyG9uNZQHf1CCOVWK3xPdbcaQuzMKiM00IeZuoZM52yD8dc7PZP+WjjCEBwFYoUQ44UQ3sA9wCeif4QQM4F/ooxARaftI4QQPvbXocB1QLYDNA0pskrquFjbpO+IKG8HIB3WoL6/LJkSgadFkKprC8v4tWBtVWWMhwgt7Vb25VawPCEcDx1DpqtOQ80Zw7/7AzYEUsp24DFgJ5ADvC6lzBJC/FwI0REF9DsgEHjjqjDReCBNCJEOvAc8KaU0DUEfSc0uxyJgaXy40VL6R+52GBED4VcHm7mWYH8v5k8Yyc4sTdcJxsxT5YuHUBjpR2eqaWi1siJB02ihXPv/1RRj1sY6cIgfQUqZAqRcte3HnV4v6+a4j4AkR2gYyuzKLid5XAgjA32MltJ3muvg7H6Y+6hybxjMisQIfrwli4KKK0wK1ywCxeKhSnNkb1EJSgYkJrma1KwyArw9WDjJtSHHDiN3O0TNVhniBqJhnKFJZ4pqGskprdM3Wqhgt3JnGDw17qDj76jtrCD+FlWv5tz7RitxOlabWhtbHBeOj6fravc7jLoS1XvDDb77piHQnA5/traGIHc7+I90WZPunhgV7Mf06GB91wnG36gSk4aAe+hk0SWqrrRqvDZmd6IYHCQBpiHQnl3ZZUyJGEaMrvVV8lOVO8OF3Zh6YkViJOlFtZRdbjZaSt/x8oVJy9RNxqZhlnQf2JlVjpeH4KY4jdfGRk6C0MlGKzENgc5camjlyNkafWcD5w4oN4YbjIg6szJR/T13aZtcdgtcKVeJSoMUKSU7s8pYMDGUIF/X1u53CE21cPZ95RZyg7Ux0xBozN7cCmxSLXBqSV4KePnDhMVGK/kEE8MCmRAaoK97KHa5SlDKecdoJU4jv55Py0gAACAASURBVOIK56sbPzba2lGwG2ztMMX49QEwDYHWpGarlpRJUcYlovQbmw1yU2DSUvDyM1rNJxBCsCIxkoNnqrnc1Ga0nL7jGwzjb1ChiVLDLOlesNPeknJ5vKaGIHcbBIRDdLLRSgDTEGhLc5uV909XsTwhQs/6KqUnoL7E7dxCHaxIjKDdJnlP1xaW8WuhphAqc41W4hRSs8uZOXY44UHGFGkbEO0tkL/brdbGTEOgKR/kV9HUZtXXLZS7XfUeiF1htJIumRE9nPBhPvoWoZuyGhCDMnqopLaJjIuX9S05ffYAtNa71SDINASakppdxjBfT+aN1ziRJuY6l7Wk7CsWi2B5QgT78ippbrMaLafvDIuE6DmQO/jWCXbZ1260DRvN3aZCfMffYLSSjzENgYZYbZI9ORXcNCUcb08N/wurCpTLwk0WyrpjRWIkja1WPiyoMlpK/4hbA6XpcLnYaCUOZWdWGbHhgUwI0yzzG+y9B1JUiK9BvQe6QsO7iMnxC5eobmjV1y2UZ2zvgd6yYMJIhvl46ptl3OF6yB08PQpqG1s5rHPI9MVjKrTXjdxCYBoCLUnNKsPbw8KNkzXty5C7HSKnqfr5boy3p4XFceHsyanAqmOPgtBJEDrlv4XNBgF7c9X/xQpd1wdyt4HF09DeA11hGgLNkFKSml3OgokjGaZjIk19ORQdcbsRUXesSIiguqGVExcuGS2lf8StgfMfQpOm+q9iV3Y5EUE+TNMxZBrsa2PXg99wo5V8AtMQaEZHIo22bqHT7tF7oLcsnhKGl4fGPQri1qjEpfxdRisZMM1tVvafrmRZfAQWHXsPVJ6G6ny3/O6bhkAzUrN0T6TZDsPHQUSi0Up6xTBfLxZMDCU1qwypY3LW6FkQGDko3EMHz1TT2GrVd33ATXoPdIVpCDQjNbucGWM0TaRpqYfCfcotpFES3IqECM5VN1JQccVoKX3HYlGL8vm7oU3DInqdSM0uI9DHkwUTNQ6ZHj0LgqOMVvIpHGIIhBA3CyHyhBAFQojHu/jcRwjxmv3zw0KImE6ffd++PU8IsdIRegYrpZebOFV8WV+30Me9B9xvRHQtOkagWruH2hpUAyBNsdkku3MquHFKmJ69B+rLVBFAN3QLgQMMgRDCA/grsApIAO4VQlzdc/AR4JKUchLwR+A39mMTUD2OE4Gbgb/Zz2fSBbs/TqTRNWIiBfxCYMx8o5X0iYggX6aPGa6vIYi5HryHqRGpppwsrqWyvkXfJLKPew8MUkMAzAUKpJSFUspW4FVg/VX7rAc22l9vApYKVSBnPfCqlLJFSnkWKLCfzzkc+D3s+onTTu9sUrPLmRAaoF8LRQBrG+TvVPVVPBzSIdWlrEiIIL2olvI6Dd0rnj4qXDEvBWwaZkkDqVnleFoEi6do3HsgZAKExRmtpEscYQiigKJO74vt27rcx97s/jIwspfHAiCEeFQIkSaESKusrOyf0ppCSHtONUTRjMtNbRw8U81yXd1C5z+E5stuuVDWGzpGort0nRXErYGGSijWs0fBruwy5k8YSbCfhiHTzXWq98CU1W67NuYIQ9DVv+zq8Iru9unNsWqjlE9LKZOllMlhYf1MpIpbCy2X4fwH/TveQPblVdBukxq7hbaDpx9MXGK0kn4xKTyQmJH++rqHOnoUaBg9dKbyCmcqG/SNFvp4bcx9c2ccYQiKgTGd3kcDJd3tI4TwBIKBml4e6zgmLFaNUDRMuU/NLic00IeZY9wrEaVXSKn+5hNvAm9/o9X0i//2KKiivlnXHgXXK4OsWRhsxyxsma6GIC8F/ENhjPO83gPFEYbgKBArhBgvhPBGLf5uvWqfrcBD9td3AHulCsreCtxjjyoaD8QCRxygqWu87CNSzX4MLe1W9udVsjwhXM9EmrJTUFfstgtlvWV5QgRtVsm+vH66Jo0mbg3UnIGq00Yr6RO7ssuZGhVE1HD3amDUK6xtcDoVptzsNr0HumLAhsDu838M2AnkAK9LKbOEED8XQqyz7/YsMFIIUQB8C3jcfmwW8DqQDbwLfEVK6dzVrLi1qiFKyQmnXsaRHDxTzZWWdr3dQsICk282WsmAmDV2BCMDvPVdJ+hYn9HIPVRZ38LxC5dYHq/pd//cB8od7eaVdh0SviGlTAFSrtr2406vm4E7uzn2CeAJR+joFZNXqoYoeSkQNctllx0Iqdnl+Ht7aJxIk6JCRgNCjVYyIDwsgqXx4ezIKKO13aZfCfCg0SqhKXc7XP9to9X0ir255UiJvusDuduVO3riTUYruSaafZMdgH8IjFuoTUy1zSbZnV3O4ilh+Hq579SyWy6dg/IM7ZLIumNFQiT1Le0cPltttJT+EbdGlUKuKzVaSa9IzSonargf8aOGGS2l70ipBpwTl7hdX+6rGXqGANSPoSIbqs8YraRH0otrqahv0dctlLdDPWsaNno1i2JD8fPyIDVLU/dQR+RKnvsHTDS2tvNBQRUrEnXty30S6i5qsTY2NA1Bx01Jgx/DruxyPCyCm3ROpAmLh5ETjVbiEHy9PLhhcii7ssv1LEIXNgVCJmoxI37/dBUt7TaN3UIpam0s1v0r5wxNQzBiHEQkaRFGmppdzvwJIQT7a5hI01gD5z8aNG6hDlYkRFJW10zGxctGS+k7Qqj/j7PvqwQ/NyY1u4xgPy/mxrhnX+seyd0OYxdCgPuv7Q1NQwBqulZ0CK64byhgYeUVCiqu6OsWOr0TpFWLqXFfWBIXjodF6O0esrWpRCc3pd1qY29uBUvjwvH00PA2VXMWKrK0GQRp+Bd2EHGrQdrg9LtGK+mWVO0TabbDsFEwaqbRShzKiABv5sSM0DeMNHoOBIS5tXso7fwlahvb9HULdbidNVkbG7qGIHIaBI9x63WCnVllJEUF65lI09YEBXvVD8Ey+L5myxMiySuv53x1g9FS+o7FQxX/O50K7S1Gq+mS1KxyvD0t3KBtX+4UCE+EkPFGK+kVg+8X2luEUC6LM3uh1f1+zOV1zZy4UMtKXYvMFe5XNfA1mRr3Fe2L0E1ZA631cO6A0Uo+hZSSXTllXDdxJAE++lWqpaEaLui1NjZ0DQGo0Wp7szIGbkaHW2hloqbrA3nbwScIYm4wWolTGBPiT1zkMH3XCSbcCF4Bbukeyiuvp6imiRW6fvfzdyq3s0ZrY0PbEIxbCL7D3TJ6KDWrjPG69h6wWVX+wKRl4OlttBqnsSIxkrTzNVRfcU/3yjXx8oNJS9X/k81mtJpPsCurHCFgabzGIdNBUTBqhtFKes3QNgQeXqr+zekdYG03Ws3HdPQe0DaRpjhN1b7XaETUH1YkRGCTsCe3wmgp/SNuDdSXul3drY/7cg/TsC93ayMU7HHr3gNdMbQNASg/XtMluHDQaCUf816u6j2grVsod5uqfR+73GglTiVxdBCjg331XSeIXaHqbrlREbrSy01kXLysb7RQ4T5ob9JqfQBMQwATl4KHj1tFD+3MKiN8mA8zojXsPQDqbxmzSNXAH8QIIVieEMGB/EqaWjVsAekfAjHXudU6gfZ9ufO2g08wjFtktJI+YRoCn0BVGTB3m1v0KGhus7Ivr5IViRF69h6oPA3VBYPeLdTBisRImttsHMh338TEaxK3FqryoKrAaCWA5n25bVbIe1fNhDVbGzMNASh/Xu0FKM80WgkH8qtoarPq7RYCbRJpBsrc8SEE+Xrq28Jyyir1nGf8rKCuuY1DhdX6uoWKjkBjlZaDINMQgP3HINwieig1q4xhvp7MG+/+9Um6JC8FRs+E4CijlbgELw8LS+LC2ZNTTrvVvaJvesXwsSq50g3cQ+/lVtBmlazQNXemY21s0jKjlfSZARkCIUSIEGKXECLf/jyii31mCCEOCiGyhBCnhBB3d/rsP0KIs0KIk/aHMfFWgeEwZp7hi2btVhu7c8pZGheuX9MTUDXui4+6fTcmR7MiMZJLjW0cO3/JaCn9I26tGs1eMTb66d1MtTY2c8ynbiPuj5TKmI6/AXyDjFbTZwZ6t3kc2COljAX22N9fTSPwoJQyEbgZeEoI0XkV9DtSyhn2x8kB6uk/catVb93aC4ZJOHruEpca2/R3CyWsu/Z+g4wbJofh7WHR1z0UtwaQ/+0dYQBNrWptbGVipJ5rY+WZcOksxN9itJJ+MVBDsB7YaH+9Edhw9Q5SytNSynz76xKgAnC/AiIdDTsMnCLvzCrDx9PCjVPc78/TK3K2QuhkVfN+CBHo48l1k0bq26MgIlG5iAz87u8/XUlTm5Wbp2o6CMp5R/Ue6LiPaMZADUGElLIUwP58zVRAIcRcwBvo3BrsCbvL6I9CCJ9rHPuoECJNCJFWWemECI2REyE8Qf2HGoCUkl3Z5VwfG4a/t6b1Vc59CPFDazbQwfKESC7UNJJbVm+0lL4jhLqBFe6DFmP078wqY7i/F/PGa9p7IHur6j0QqOcgrkdDIITYLYTI7OKxvi8XEkKMAl4APiul7FhV+z4QB8wBQoDvdXe8lPJpKWWylDI5LMxJf+z4daqRSr3rp/iZF+u4WNukb5G5vO2q98AQcwt1sDwhAiGUn1tL4taAtUVlxbqY1na1NrY8PkLP3gNV+VCZo61bCHphCKSUy6SUU7t4bAHK7Tf4jht9l6tNQoggYDvwQynloU7nLpWKFuA5YK4j/lH9JmEdIA1ZNE7NLsMiYGm8poYg553/RqAMQcKG+TA3JoQdmXo0hf8UY+aDX4ghiZUfnamivrldY7fQVvU8mA1BD2wFHrK/fgjYcvUOQghv4C3geSnlG1d91mFEBGp9wdhA/vAEGDnpv/+xLmRnVhlzx4cQEqBXIgqgWh6eeU/NqDSqr+JoVk2N5HS56iqnHR6e9rpb74K1zaWX3plVZl9nCXXpdR1G9laIStY6ZHqghuBJYLkQIh9Ybn+PECJZCPGMfZ+7gBuAh7sIE31JCJEBZAChwC8HqGdgCKFuZmcPqH67LuJsVQOny6/oGy10eqdqfZjQJ2/hoOPmqaMAeFfXWUH8WmXUz77vsktabZLUrHJuigvH18vDZdd1GLUXoPSk1rMBGKAhkFJWSymXSilj7c819u1pUsrP21+/KKX06hQi+nGYqJRyiZQyye5qekBKafxQKmGd8nW7MIKiw6+sbf31nK2qJWVUstFKDCUy2JdZY4eTkqHpOsHEJeAdCNlvu+ySR8/VUN3Qyipt3UL24BLN18Y0XJlxMqNmKF+3C91DKRmlTI/WtCVlawPk71ZRJ4OwJWVfWZ00iuzSOj1bWHr5weSVkLPNZWXZ3820h0zr2pIyeytEJEHIBKOVDAjzl3s1He6hM++pabKTuVDdSMbFy6yZNsrp13IKBbtV2V3NR0SOomPBc4eu0UMJG6CpxiUtLG02yc6sMm6YHKZnS8r6Mig6rL1bCExD0DUJ65XPO+9dp18qxe5PXjVVU0OQ846KNhm70GglbkH0CH+mRQfrawhil6sWli5wD526eJnSy836uoVytwFyUAyCTEPQFVHJyuftAvfQ9lPKLTQmxN/p13I47S1qoThujYo6MQGUUU8vqqX4UqPRUvrOx+6hd5zuHtqRWYqnRbA0TtOQ6eytMDIWwuKMVjJgTEPQFRaLmu4V7IYW561fd7iFVidpOhso3ActdUM+WuhqOka42iaXJW6Axmo4/4HTLiGl5N3MMhZMHEmwv5fTruM0Gmvg3AfqPjEIQqZNQ9AdCeuhvRkKdjntEh1uIW0NQfYW8AmC8TcarcStiAkNIH5UkL6GYNJy8PJX/79OIuPiZc5XN3LLtNFOu4ZTyR1cmfSmIeiOsQsgIMypPwbt3UI521S0kGbdmFzB6qmRpJ2/RHlds9FS+o63v+pnnPOO6rrlBLadKsXLQ+ibO5O1GUbEqCjDQYBpCLrD4qF836dTodXxvl7t3UIFe6DlMky93WglbsmqJHWD25ml6awgcQM0VML5Dx1+aptNsv1UKdfHhunpFrpSCYX71Xd/ELiFwDQE1ybxVmhrgPydDj+19m6hzDdVtNAE0y3UFZPChxEbHsj2U5pmGceuAE8/yHJ89NCJoktcrG3ilumafvdztii30CAaBJmG4FrEXA+BEZCxyeGn1tot1NqompgkrAMPDUd0LmLNtFEcOVejqXsoACY7xz30Tnop3p4WlulaYDFzM4ROUbXJBgmmIbgWFg81K8jf5dDkMu3dQvk71UxpEI2InMHaaaOREn1nBQkboKFClWZ3EFabJCWjlJumhDHMV8NBRF2J+nsMIrcQmIagZ6beoeq0O7D20KBwCwVGwLjrjFbi1kwKDyRhVBBb00uMltI/OtxDDkwuO3K2hor6Fm6Zrmm0UNbbgISptxmtxKGYhqAnopNV7SEHuoe0dgs116kF9MRb1YzJ5JqsmzGak0W1FNVomFzmE6gyjbO3OCy5bNupEvy8PFgSd81mhu5L5puq50ZorNFKHIppCHpCCDUNLNwHDVUDPt25qga93UJ5KWqGZLqFesUa+//zO6c0nRUk3amih87uH/Cp2q02dmSWsTQ+XM92rJfOwcW0QfndNw1Bb5h6h4oSyHprwKfqcBNoOzXOfBOCx0D0HKOVaMGYEH9mjR3OO+marhPErlBJgw6YEX90ppqahlbW6ppElrlZPSfeaqwOJzAgQyCECBFC7BJC5NufR3Szn7VTU5qtnbaPF0Icth//mr2bmfsRkajqiWS+OaDTSCl5++RF5o4PYbSOJacba+DMXvVDGEQLZc7mlumjySmto6BCw8b2Xr6qGm/OO9DWNKBTbTtVQqCPJ4unaFpyOnOzGgCNGGe0Eocz0BnB48AeKWUssMf+viuaOjWl6ZyT/Rvgj/bjLwGPDFCPcxBCzQouHITaon6fJqukjsLKBjbM0LSlXc47YGsflFNjZ7Jm2igsArbqOitIugNa6yE/td+naG238W5mGSsSIvTsRFaZB+UZg/a7P1BDsB7YaH+9EdV3uFfY+xQvATrmnH063uV0RAlkbe73Kd4+cREvD8HqJE3T6jPfhJCJMGq60Uq0InyYL/MnjOSd9BKklEbL6Tvjb4CAcMh4o+d9u+FAfiV1ze2s1TWJLHMzIFRI7SBkoIYgQkpZCmB/7i4UwFcIkSaEOCSE6PhLjgRqpZQd4QjFQLdDZSHEo/ZzpFVWVg5Qdj8YORFGz+q3r9Rqk2xNL2HxlHCG+7unB+ya1JerZiVTbzPdQv1gw4wozlY1cLKo1mgpfcfioUbCp1OhqX/63zpxkRH+XiyapKFbSEo1AIxZBEGaGrIe6NEQCCF2CyEyu3j0pfbwWCllMnAf8JQQYiLQ1d2k2+GSlPJpKWWylDI5LMygL1PSHVB2Cqry+3zo4cJqKupbWD9D14WyTSBtykVm0mdWJUXi42nhrRMXjZbSP5LutOfTbOvzoXXNbezKLueW6aPx9tQwPqXkBFSdHrRuIeiFIZBSLrM3l7/6sQUoF0KMArA/V3RzjhL7cyGwD5gJVAHDhRAdcWTRgHvH2CXeCoh+zQrePnmRQB9PfdPq01+B0TMhXP8mHEYwzNeL5QkRbE0vobXdZrScvhM1C0aMh1Ov9/nQHRmltLTbuHWmpmtj6a+Ch8+gjBbqYKDmeSvwkP31Q8CnajYLIUYIIXzsr0OB64BsqZyl7wF3XOt4tyJoNIy/Hk69qqaLvaS5zcqOzDJWJkbquVBWlgllGTD9PqOVaM3ts6KpbWxjX16X4yX3RgiYdjecfR8uF/fp0M3HLzIhNIAZY4Y7SZwTaW9Vs+G41eCnof5eMlBD8CSwXAiRDyy3v0cIkSyEeMa+TzyQJoRIR934n5RSZts/+x7wLSFEAWrN4NkB6nE+0+9TiSUXDvb6kN055dQ3t7NhpqZuoVOvgsVzUE+NXcH1saGEBnrr6x6afg8g1Qi5lxTVNHL4bA23zoxC6Li2VLBbdWubfq/RSpzKgAyBlLJaSrlUShlrf66xb0+TUn7e/vojKWWSlHK6/fnZTscXSinnSiknSSnvlFK2DOyf4wIS1oF3IJx8qdeHvJFWzOhgXxZODHWiMCdhbVfugNiVEDDSaDVa4+lh4Zbpo9mTU8Hlxjaj5fSdkPGqvlT6K72eEW85qYzeBm3dQi+rBlUTlxitxKlouHJjMN4BKoQs621obehx97LLzRzIr+T22dF4WDQcERXugyvl9tGgyUC5bWY0rVYb2zLcezmsW2bcB9UFUHy0x12llGw+oRIotayr1VgDee9C0l2Dvty6aQj6w4z7oPWKatXYA5tPFGOTyj+sJemvgO9wmLzSaCWDgqlRQcSGB7L5uKbuoYT1qp9xL2bE6cWXKaxs4PZZms4GsjaDrW1IDIJMQ9Afxi5Q/Up7+DFIKdmUVszcmBBiQgNco82RNNepcMGpt4Onj9FqBgVCCG6fHc2x85c4U3nFaDl9x2eYKjmR+VaPJSc2HSvCx9PCKl0LLJ58BSKmwqhpRitxOqYh6A8Wi1o8Ovv+NUtOHL9QS2FVA3fM1nQ2kPkmtDerGZCJw7htVhSeFsHrR/tfrsRQZtyn+lVfo0dHY2s7W06UsCZpFEE6NqCpzFOVRofAbABMQ9B/ehFBselYEX5eHqyepumI6PhGCE+EqNlGKxlUhA/zZWl8OG8eL9YzpyDmelWB9hoz4pSMMupb2rl7zhgXCnMgxzaCxWvQRwt1YBqC/jIiRv0gTrwAtk//mJtarWxLL2VVUiSBPhrWXi9NVxmVsx8yS0o4gXvmjKXqSit7c8uNltJ3LBY1KzjzHlw63+Uurx65wITQAOaOD3GxOAfQ1qyiheLXQoCGkX79wDQEAyH5s1B7Hgr3fuqjbadKqG9p565kjUdEnr4w7S6jlQxKbpgcRmSQL68c0dQ9NOtBNUA4vvFTHxVU1JN2/hJ3zxmjZ+5A7jZougSzHup530GCaQgGQtwt4B8Kac996qMXD19gUngg83QcEbU2qNyBhA3g12WLCZMB4mER3JUczfv5lVysHVidf0MIjlZNa068CNZP5kS8eqQIT4taFNeSY/9RM/7xNxqtxGWYhmAgeHrDzPshbwfU/TcuPPPiZdKLarl/3lg9R0RZb6n687OHzojICO60zxbfSNN0VjD7syrHJC/l400t7VY2n7jI8oQIQgM1jDSrKlBVdmc9qFxgQ4Sh8y91FrMfVm0sj7/w8aaXDp/H18vCbbrmDhzbCKGTVZisidMYE+LPokmhvH60CKtNwz4FscshKBrS/v3xpl3Z5dQ0tHLP3LEGChsAxzeqciozHjBaiUsxDcFACZkAE25SXyBrO3XNbbx9ooR100cT7Kdh2Fx5FhQfUf5RHWczmnH/vLGUXG5md46Oi8YeatZYuA+qzwDwwsHzRA33Y9EkDRdZ21vh5Msw+WYYpmmV4H5iGgJHkPw5qLsIBbt4+8RFmtqsPDBf076mR55Wi8RDJGzOaJbFRzA62JeNH50zWkr/mPkZEB5wfCPZJXUcPlvDQwvH6VlOJestaKxSQSBDDNMQOIIpqyAwAnn0WV48dJ6kqGCmRWtYsraxBtJfU5FCZoE5l+DpYeEzC2L46Ew1eWUaNrcPGqW+/yde5KUP8vDz8uDuZA3dQlLC4b8rl+jEpUarcTmmIXAEHl6Q/AiiYBfWijwemK/hDwGUe6u9CeZ9yWglQ4p75ozBx9PCxoPnjJbSP+Z9UZVqznidDTOjCPbX0CVadETlzcz74pB0iZqGwFEkf4424cWXfVNZN13DIlvWdjjyjEqSi0g0Ws2QYkSAN+tnjOat4xf1LE8dcz1VAZN5SGzn4QWaukQP/x18g4esS9Q0BA7iTJMfm9oWsUG8j1+bhg3Kc9+BumKY/2WjlQxJHloYQ1Obldc1DCVtt0n+1rKSyZaLTGk4YrScvnO5GLK3qpBRbw2LQzqAARkCIUSIEGKXECLf/vyp7CMhxE1CiJOdHs1CiA32z/4jhDjb6bMZA9FjJM9+cJYXWI2nrQWO/bvnA9yNQ/+A4eNUxISJy0kcHczcmBCeP3ROu1DSXdnlvHBlDi2+YXDwb0bL6TtH/gVImPuo0UoMY6AzgseBPVLKWGCP/f0nkFK+J6WcIaWcASwBGoHUTrt8p+NzKeXJAeoxhOorLbx5rJhpM+fDpGXqi9Xu/s3WPubicSg6pPyjFg17Kg8SPrcohqKaJlIySo2W0muklPzz/UIiQ4bhNf9ROLMHKnKMltV7WhvV2ljcWhiu6dqeAxioIVgPdBQb2Qhs6GH/O4AdUsrGAV7XrXjh0Hla2m18/vrxsOArKtsy802jZfWeD/4IPkEwc2gl0bgbKxIimRgWwN/2nUH2shWk0RwsrOZkUS1fvGEiljmPqNDjQxrNCo5vVHWFFnzFaCWGMlBDECGlLAWwP4f3sP89wCtXbXtCCHFKCPFHIUS3OelCiEeFEGlCiLTKysqBqXYgzW1WXjh4npumhDEpfJhKLgtPgI/+3GVVUrejIgdytqrZgG+w0WqGNBaL4MuLJ5FTWse+PPf5jl+Lv75XQNgwH9VzI2CkWmxNfw3qNJjVtLfAh/8H4xbB2PlGqzGUHg2BEGK3ECKzi8f6vlxICDEKSAJ2dtr8fSAOmAOEAN/r7ngp5dNSymQpZXJYWFhfLu1U3kgrorqhlS9cP0FtEAIWfQsqstUCrLtz4A/gFQDzzEVid2D9jNFEDffjL+8VuP2s4MSFS3xYUM0Xrh+Pr5fdpXjd18HWrm6w7s7Jl6C+FG74/4xWYjg9GgIp5TIp5dQuHluAcvsNvuNGX3GNU90FvCWl/Dg+TkpZKhUtwHPA3IH9c1xLc5uVv7xXwJyYESyY2CkBa+ptMDIW9v/WvWcF1Wcgc5PKpDQTyNwCLw8Lj94wgWPnL3HkbI3Rcq7J3/adIdjPi/vmdQoZDRmvZgXHnoP6MuPE9YS1TblEo5JhwmKj1RjOQF1DW4GOEpUPAVuuse+9XOUW6mREBGp94Qh5qQAAEMdJREFUIXOAelzKK0cuUF7XwjeXT/5klVGLB9z4XSjPVLXN3ZUPn1JdmBZ+1WglJp24e84YQgO9+du+M0ZL6Za8snp2ZZfz8MKYTzdeuuHb6kbrzrOCjDeg9gLc8J0hmUB2NQM1BE8Cy4UQ+cBy+3uEEMlCiGc6dhJCxABjgP1XHf+SECIDyABCgV8OUI/LaG6z8rd9Z5g3PoSFE7sosDX1dhg5Cfb/xj1nBbVFqjn3rAdhWKTRakw64evlwecWjWf/6UqOX7hktJwu+fPefPy9PXh4YcynPwyZoFq5pv0b6t2wmJ7NCgd+DxFJMHml0WrcggEZAilltZRyqZQy1v5cY9+eJqX8fKf9zkkpo6SUtquOXyKlTLK7mh6QUl4ZiB5X8uKh81TWq9lAl1g84Ab7rCCv+ybfhrH/N+r5uq8bq8OkSx5aEENooA9PpuS63VpBelEt206V8rnrxjMiwLvrna5341nBqdehukCtDZizAcDMLO4Xja3t/GP/Ga6bNJL5E67hW++YFbz3KzUKcRfKMlVnqXlfhOGattIc5AT4ePKNZbEcOVfDnpxrLb25FiklT+7IJSTAmy/eOKH7HUdOhGl3Q9qzcPmi6wT2RGsj7P0FjJ4F8euMVuM2mIagHzx74CxVV1r55rJuZgMdeHjCkh+qCKIuersaxq4fg2+QGrWZuC13zxnDhNAAfvNurttkG+87XcnBwmq+vjSWYb49FJdb/Liq6rn7py7R1isO/VWVjF/xyyHVgawnzL9EHympbeKv+wpYnRRJckwv+hEnbICxC2HvL6HJDWoQFexR2Z83fBf8NeynPITw8rDwnZVTyK+4wpvHio2Wg9UmeTIll3Ej/bm3Nx3IRoyDhY9BxutQdNT5AnviSgV88JTKIo65zmg1boVpCPrIr1JykBJ+sDq+dwcIAaueVLX+3/+dc8X1hM2qZgPDx8HcLxirxaRX3Dw1khljhvOHXadpajXWvbj5eDF55fV8d2Uc3p69vHUs+hYERsK7jxsfNLHv19DeDMt+ZqwON8Q0BH3g4Jlqtp0q5cuLJxI9wr/3B46arso3HP6Hao5tFCdfUovXy34Knho2Fh+CCCH43zXxlNU185f38g3TUdfcxm935jFjzHBWJ/UhyswnEJb9BC6mqZwVo6jIUb24kx+B0EnG6XBTTEPQS9qtNn72ThZRw/340o0T+36CpT8GTz/Y+QPlN3U1DVVqNjB2ASTe6vrrm/SbOTEh3DYriqffL6SgwpjAut/vzKP6Sgu/WD/1kzkzvWHaPTB6Juz6CbQY0IXNZoNt31TrYjd2W7xgSGMagl7yn4/OkVtWz4/Wxv83nb4vBIbD4u9B/k7VG9XVpP4QWq7A2qfMkDkN+cHqePy8PPjBWxnYXLxwfLKolucPnecz88eRFN2PelQWC6z6nSrnsOfnjhfYEyeehwsH1QKxmUHfJaYh6AWFlVf4f6l5LIkLZ2XiAJKv5n1Zha2lfAcaqh0nsCfy3oX0V1TOQHic665r4jBCA3344doEjpyt4fmD51x23eY2K99+/SSRQb58e+WU/p9ozBzVAvXI03D+I8cJ7InaIkj9kSosN+N+111XM0xD0APtVhvf2XQKbw8Lv74tqe/T4s54eML6v0LzZXjna65xETVUw9avQsRUVfbCRFvunB3N4ilh/ObdPAorXeMi+n1qHmcqG/jN7dMI6ilctCeW/FAFKrz1JWiuc4zAa2GzwdtfBmmD9X8xZ8LXwDQEPfDU7nyOnb/ELzZMJSLId+AnjEhQi2e521RhLmdis8GWr6h667f+01wg1hwhBE/eNg1fLwuPvXyC5jbnRhG9l1fBvw6c5f55Y7lhsgMq/voEwm3/Uq0ht7sgh+WjP8G5A7DyV6oYnkm3mIbgGuw/Xclf9xVwV3I062c4sCH9/K/AxKXw7veh5ITjzns1H/0fnN6hfgiRU513HROXERnsy+/vmk52aR2/2JbttOuU1Dbx7dfTiYscxo/WJjjuxGPnqUSzjNfh6LOOO+/VnPsA9vxM5fHMetB51xkkmIagGwoqrvDYy8eZEjGMn61z8E3UYoHbnoaAcHj1fpXo4mhOp6qFucTbzJyBQcaSuAi+eOMEXjp8gRcOnnP4+Rta2vn8xjRa22389f5Z/QuOuBbXfxsmLYcd33XOesGlc/DGw6r43bo/my6hXmAagi6oqG/mkY1H8fG08MxDyfh5O6GPb0Ao3POSSjR76Q7H+kxL09UPITLJ/CEMUr67Mo6lceH89J1s9uY6rsJnm9XG1189QW5ZHX+5byYTwwIddu6PsXjA7c+o9YJX7nVsj+PGGnjpTlXw7p5XVMioSY+YhuAqahpaeeCZw1TWt/D0g8l9SxzrK6OmwV3PQ3mW+kG0OGABsCwTnt+gykfc97ryy5oMOjwsgqfumUHCqCC+9OJxPiqoGvA5rTbJN187ye6cCn62fiqLp/TUeXYA+A2HB95U61Yv3KaaJA2UpkvwwgY1I7jnZQjroRaYyceYhqATRTWN3PmPjzhf3cgzDyYza+wI51908gq1kHvhIDy/Xo1o+suFw7DxFvDyg4e2mn0GBjnDfL3Y+Lm5xIz057P/Ocq7mf3vCNbUauVLLx5j26lSfrA6js/M///bO9PYKq4rjv9O/LDZF9sN++YS1gAiJYBRmiYkBBdSaAQpNKUl0IgqhUpVS5sQ+gFaVUqjKpHyoU0c0YREpCyBqKikQYBaiigyAcISgwAHm80EMGY14PX0wx3ah3nm2Z63+b3zk57ezJ17Z87/zbw5M3Pu3NM3fCO/ZPaH2euhthL+kgdfHWz+uq6WumP//GGYudLGEmoi5gg8dhSV8cyf/sOFa5W8P28M4weESDYTLYbPgJkfuD9C/mNQuq9p7VVh7wew4ml3pfX8393zUSPpyWyXzqr5uQzt0ZEXV+7hza3HmjxS6cmLN5iZv5Mth8+xbOow5j/ajDfnm0u3B2HuPyCtFSyfBAebMQzFyQJ4ZwKUF7vHQQOfirydSY4vRyAiz4pIoYjUicjoe9TLE5EjIlIkIi8HlfcXkQIROSYiq0WkgSwX0ePyjSqWbihk9vICOrUJsO7F8Yy9V46BaDF4Cjy/0SX+Xj7RjVZadSN8u/JiF3DesBB6j4UXtpoTSDEy26Xz4QvjmDqyB69vPsqs/J0Ull4J2666to73dhQz+c3tlJRV8PbsbzAnVMaxaPO1Qe647TYc1v0Y1s6Fq2fDt7t1xb0s9m6ecyTzNsEDT0bf3iRE/GQ/EpEhQB3wNrBIVXeHqJMGHMWlsjwNfAZ8X1UPicgaYL2qrhKRt4D9qvrncNsdPXq07t5916YaTU1tHV+UXmXDvlLW7jlFRWUNPxjbl8WTB9M2PRB+BdGk4iJsWgwHVkPbbBg9D4Z8x70Qdnv89KobcKrAZVo6uAbu8/IejFtgY6ynMKrKur1n+P3GQ1y+Wc2TQ7ryvdG9GZeT+b/cAarK8bIKNh86x8qCE5wqv8kjA7J5dfrw6MbDGkNtNWx/3aWRBBg5E4Y/C73HQcC7RqytcQMnFq53g8jdugyjfui6SFtgOCwiskdV77po9+UIglb+Lxp2BLnAUlWd5M0v9ha9ClwAuqlqTf1696K5juCVjw+y7cgFyq5XUllTR6s0YdKwbix4fABDuifYQXRip0suf/RTNx9oDe27ujuGa2fd25IZHV0WqG/+Ejp2j6+9RsJw5WY17/z7OH/ddZKLFVWIuCEq2rRK41JFFdcqawAY0y+Tn3wrhwmD7/f3xnykKS92L4Pt+9ANGy1p0LGHW1ZxwZXdF4CBee5t+e4j42tvC6IhRxCLy9+ewKmg+dPAWCALuKyqNUHlDb61JSLzgfkAffo0IilGKEM6t2Fs/0yyO2QwolcncnOyyGqfoG/b9s11n2vnoGiLy3J2/Tykpbs/Rc+HIOcxFxg2jCA6tWnFokmD+NkTA9hTcoldJeV8deUWN6tr6dI2nUHdOpCbk0W/7HbxNjU0mf3h6Tdg4u+geBuc2eveRhaBtlnuxP/1Ca4LthERwjoCEdkChOp+skRV/9aIbYS61NB7lIdEVfOBfHB3BI3Y7l0seLwFjkPeoSuMssGyjKaTEUhj/IDs2HZ8iCQZ7V3sbPCUeFuS9IR1BKrqN/pyGgjOkN4LKAXKgM4iEvDuCm6XG4ZhGDEkFpHFz4AHvB5C6cAsYIO64MQ/gRlevTlAY+4wDMMwjAjit/voMyJyGsgFNorIJq+8h4h8AuBd7S8ENgGHgTWqWuit4iXgFyJShIsZRHEUKsMwDCMUEek1FGv8dh81DMNIRRrqNWSdzg3DMFIccwSGYRgpjjkCwzCMFMccgWEYRorTIoPFInIBONHM5tm4dxhSCdOcGpjm5Mev3r6qelcC6hbpCPwgIrtDRc2TGdOcGpjm5Cdaeu3RkGEYRopjjsAwDCPFSUVHkB9vA+KAaU4NTHPyExW9KRcjMAzDMO4kFe8IDMMwjCDMERiGYaQ4SeUIRCRPRI6ISJGIvBxieYaIrPaWF4hIv6Bli73yIyISNl1mItBcvSIyUUT2iMhB73tCrG1vLn72sbe8j4hcF5FFsbLZLz6P6xEislNECr393TqWtjcXH8d2KxFZ4Wk9HJQaN+FphOZHRWSviNSIyIx6y+aIyDHvM6fJG1fVpPgAacCXQA6QDuwHhtar81PgLW96FrDamx7q1c8A+nvrSYu3pijqHQX08KYfBM7EW0+0NQctXwesxeXYjrumKO/nAHAAGOnNZyX6cR0Bzc8Bq7zptkAJ0C/emiKkuR8wAngfmBFUngkc9767eNNdmrL9ZLojGAMUqepxVa0CVgHT6tWZBqzwpj8CnhCXtXsa7uCpVNVioMhbXyLTbL2q+rmq3s4GVwi0FpEETd58B372MSLyXdyfpJCWgx/NTwEHVHU/gKpeVNXaGNntBz+aFWgnIgGgDVAFXI2N2b4Iq1lVS1T1AFBXr+0kYLOqlqvqJWAzkNeUjSeTI+gJnAqaP+2VhayjLmHOFdxVUmPaJhp+9AYzHfhcVSujZGckabZmEWmHS4S0LAZ2RhI/+3kgoCKyyXuk8OsY2BsJ/Gj+CKgAzgIngT+qanm0DY4Afs5Bvs9fYXMWtyAkRFn9vrEN1WlM20TDj163UGQY8AfclWNLwI/mZcAbqnrdu0FoKfjRHAAeAR4GbgBbvcQkWyNrYsTxo3kMUAv0wD0m2S4iW1T1eGRNjDh+zkG+z1/JdEdwGugdNN8LKG2ojnfr2Akob2TbRMOPXkSkF/Ax8CNV/TLq1kYGP5rHAq+JSAnwc+AVEVkYbYMjgN/jepuqlqnqDeAT4KGoW+wfP5qfAz5V1WpVPQ/sAFrCWER+zkH+z1/xDpJEMNgSwD3/7c//gy3D6tVZwJ0BpjXe9DDuDBYfJ8GDaj71dvbqT4+3jlhprldnKS0nWOxnP3cB9uKCpgFgCzAl3pqirPkl4F3cVXI74BAwIt6aIqE5qO573B0sLvb2dxdvOrNJ24/3DxDhH3MycBQXfV/ilf0WmOpNt8b1GCkCdgE5QW2XeO2OAN+Ot5Zo6gV+g3uOui/oc3+89UR7Hweto8U4Ar+agdm44PgXwGvx1hJtzUB7r7zQcwK/ireWCGp+GHf1XwFcBAqD2s7zfosiYG5Tt21DTBiGYaQ4yRQjMAzDMJqBOQLDMIwUxxyBYRhGimOOwDAMI8UxR2AYhpHimCMwDMNIccwRGIZhpDj/Ba3SaA0burncAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "timing = np.linspace( 0, 5/signal.freqs, int(signal.freqs*4))\n",
    "raw = mysignal(timing)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(timing, raw.real)\n",
    "ax.plot(timing, raw.imag)\n",
    "ax.plot(timing, np.abs(raw))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Correlate signal with a sine wave\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ u[n] = \\exp(2 \\pi i f \\frac{n}{f_s} + i\\phi_t)$$\n",
    "sampled with $f_s$, such that $N = f_s T$\n",
    "\n",
    "$$ v(t; f, \\chi) = \\exp(2 \\pi i f t + i\\chi ) $$\n",
    "\n",
    "---\n",
    "\n",
    "$$\\begin{align}\n",
    "Corr(\\eta; u[n], v(t)) \n",
    "    &= \\sum_{n=0}^{N-1} u[n] \\; v^*(t - \\eta; f, \\chi) \\quad \\quad \\quad \\quad \\quad \\quad \\quad \\quad \\quad \\\\\n",
    "    &= \\sum_{n=0}^{N-1} u[n] \\; \\exp(-2 \\pi i f ( t - \\eta ) -i\\chi ) \\\\\n",
    "    &= \\sum_{n=0}^{N-1} u[n] \\; \\exp(-2 \\pi i f ( \\frac{n}{f_s} - \\eta ) -i\\chi ) \\\\\n",
    "    &= \\sum_{n=0}^{N-1} \\exp(2 \\pi i f \\frac{n}{f_s} + i\\phi_t) \\; \\exp(-2 \\pi i f ( \\frac{n}{f_s} - \\eta ) -i\\chi ) \\\\\n",
    "    &= \\sum_{n=0}^{N-1} \\exp(2 \\pi i f \\frac{n}{f_s} + i\\phi_t -2 \\pi i f ( \\frac{n}{f_s} - \\eta ) -i\\chi ) \\\\\n",
    "    &= \\sum_{n=0}^{N-1} \\exp(2 \\pi i f ( \\frac{n}{f_s} - \\frac{n}{f_s} + \\eta ) + i\\phi_t -i\\chi ) \\\\\n",
    "    &= \\sum_{n=0}^{N-1} \\exp(2 \\pi i f \\eta + i\\phi_t -i\\chi ) \\\\\n",
    "    &= \\exp(2 \\pi i f \\eta + i\\phi_t -i\\chi ) \\sum_{n=0}^{N-1} 1 \\\\\n",
    "    &= (N-1) \\exp(2 \\pi i f \\eta + i\\phi_t -i\\chi ) \\\\\n",
    "\\end{align}$$\n",
    "\n",
    "---\n",
    "Max correlation if $ 2 \\pi i f \\eta + i\\phi_t -i\\chi = 0 + 2 \\pi i k$ :\n",
    "\n",
    "$$ \\eta = \\frac{ \\chi - \\phi_t }{ 2 \\pi f } $$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "def correlate_function( eta, timing, raw, func, func_params=[] ):\n",
    "    N = len(raw)\n",
    "    \n",
    "    \n",
    "    np.sum( raw * func(timing - eta) )\n",
    "    \n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "$$\\begin{align}\n",
    "Corr(\\eta; u[n], v(t)) \n",
    "    &= \\sum_{n=0}^{N-1} u[n] \\; v^*(t - \\eta; f', \\chi) \\quad \\quad \\quad \\quad \\quad \\quad \\quad \\quad \\quad \\\\\n",
    "    &= \\sum_{n=0}^{N-1} u[n] \\; \\exp(-2 \\pi i f' ( t - \\eta ) -i\\chi ) \\\\\n",
    "    &= \\sum_{n=0}^{N-1} u[n] \\; \\exp(-2 \\pi i f' ( \\frac{n}{f_s} - \\eta ) -i\\chi ) \\\\\n",
    "    &= \\sum_{n=0}^{N-1} \\exp(2 \\pi i f \\frac{n}{f_s} + i\\phi_t) \\; \\exp(-2 \\pi i f' ( \\frac{n}{f_s} - \\eta ) -i\\chi ) \\\\\n",
    "    &= \\sum_{n=0}^{N-1} \\exp(2 \\pi i f \\frac{n}{f_s} + i\\phi_t -2 \\pi i f' ( \\frac{n}{f_s} - \\eta ) -i\\chi ) \\\\\n",
    "    &= \\sum_{n=0}^{N-1} \\exp(2 \\pi i [\\frac{n}{f_s}(f-f') + f' \\eta ] + i\\phi_t -i\\chi ) \\\\\n",
    "    &= \\exp{(i (\\phi_t - \\chi))} \\; \\exp{(2 \\pi i f' \\eta)} \\; \\sum_{n=0}^{N-1} \\exp(2 \\pi i \\frac{(f-f')}{f_s} n) \\\\\n",
    "\\end{align}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "Write $\\alpha = \\frac{(f-f')}{f_s}$,\n",
    "$$\\begin{align}\n",
    "0\n",
    "    &= \\frac{\\partial}{\\partial \\alpha} \\sum_{n=0}^{N-1} \\exp(2 \\pi i \\alpha n)\n",
    "    \\\\\n",
    "    &= \\sum_{n=0}^{N-1} \\frac{\\partial}{\\partial \\alpha} \\exp(2 \\pi i \\alpha n)\n",
    "    \\\\\n",
    "    &= \\sum_{n=1}^{N-1} \\frac{1}{2 \\pi i n} \\exp(2 \\pi i \\alpha n)\n",
    "    \\\\\n",
    "    &= \\frac{1}{2 \\pi } \\sum_{n=1}^{N-1} \\frac{ \\cos(2 \\pi \\alpha n) - i \\sin(2 \\pi \\alpha n)}{ni}\n",
    "    \\\\\n",
    "    &= \\frac{-1}{2 \\pi } \\sum_{n=1}^{N-1} \\frac{ i \\cos(2 \\pi \\alpha n) + \\sin(2 \\pi \\alpha n)}{n}\n",
    "\\end{align}$$\n",
    "\n",
    "Actually, we are looking at the amplitude (or power) so we take the absolute value:\n",
    "$$\\begin{align}\n",
    "    \\\\\n",
    "    &= \\frac{1}{2 \\pi } \\sum_{n=1}^{N-1} {\\frac{ \\cos(2 \\pi \\alpha n) + \\sin(2 \\pi \\alpha n)}{n}}\n",
    "    \\\\\n",
    "    &= \\frac{1}{\\sqrt{2} \\pi } \\sum_{n=1}^{N-1} {\\frac{ \\sin( 2 \\pi \\alpha n + \\frac{\\pi}{4}) }{n}}\n",
    "    \\\\\n",
    "    &= \\sum_{n=1}^{N-1} \\frac{1}{\\sqrt{2} \\pi } {\\frac{ \\sin( 2 \\pi \\alpha n + \\frac{\\pi}{4}) }{n}}\n",
    "    \\\\\n",
    "    &= \\sum_{n=1}^{N-1} \\sqrt{2} \\alpha {\\frac{ \\sin( 2 \\pi \\alpha n + \\frac{\\pi}{4}) }{2 \\pi \\alpha n}}\n",
    "\\end{align}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "We get a $\\sin(2\\pi\\alpha n + \\frac{\\pi}{4}) $ dependency, \\\n",
    "the amplitude is zero when\n",
    "$$\n",
    "    2\\pi\\alpha n + \\frac{\\pi}{4} = \\pi k \\; \\text{with k integer}\n",
    "$$\n",
    "or\n",
    "$$\n",
    "    \\alpha n = \\frac{4k-1}{8} = \\frac{k}{2} - \\frac{1}{8} \\;\\text{with k,n integer}\n",
    "$$\n",
    "\n",
    "Likewise, the amplitude is maximal when\n",
    "$$\n",
    "    2\\pi\\alpha n + \\frac{\\pi}{4} = \\pi (k+\\frac{1}{2}) \\; \\text{with k integer}\n",
    "$$\n",
    "or\n",
    "$$\n",
    "    \\alpha n = \\frac{4k+1}{8} = \\frac{k}{2} + \\frac{1}{8} \\; \\text{with k,n integer}\n",
    "$$\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\begin{align}\n",
    "Corr(\\eta; u[n], v(t)) \n",
    "&= \\exp{(i (\\phi_t - \\chi))} \\; \\exp{(2 \\pi i f' \\eta)} \\; \\sum_{n=0}^{N-1} \\exp(2 \\pi i \\frac{n}{f_s}(f-f'))\n",
    "\\\\\n",
    "&= \\exp{(\\pi i ( 2f' \\eta + \\phi_t - \\chi))} \\; \\sum_{n=0}^{N-1} \\exp(2 \\pi i \\frac{n}{f_s}(f-f'))\n",
    "\\end{align}$$\n",
    "maximal with\n",
    "$$\n",
    "n = \\frac{f_s}{2(f - f')} {(k + \\frac{1}{4})} \\text{, and }\n",
    "\\eta = \\frac{ \\chi - \\phi_t }{ 2 \\pi f' }\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "almost_sinc() missing 1 required positional argument: 'N'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-4-a147ffca3697>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     33\u001b[0m \u001b[0;31m#ax.plot(alphas, np.abs(misfit)**2, label=\"Power\")\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     34\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 35\u001b[0;31m     \u001b[0max\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0malphas\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mN\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0malmost_sinc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msinc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpi\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;36m1.5\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0malphas\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mN\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mabs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msinc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpi\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;36m1.5\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0malphas\u001b[0m \u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m \u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     36\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     37\u001b[0m     \u001b[0mn\u001b[0m\u001b[0;34m=\u001b[0m \u001b[0;36m10\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: almost_sinc() missing 1 required positional argument: 'N'"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAs8AAAEGCAYAAACafXhWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeXycVb348c+ZfTLZl6Zt0jbdtyzdoJRFNgFZvQiugFdRREQBxQWvV7j4cwEuqChwBUWKgqCyyKIsZau00JY2TZOm+5Ku2ZNmm33m/P6YTNo0yzzJTJvt+3698qrO82Seb0r6PN8553u+R2mtEUIIIYQQQsRmGuoAhBBCCCGEGCkkeRZCCCGEEMIgSZ6FEEIIIYQwSJJnIYQQQgghDJLkWQghhBBCCIMsQx3AQGRnZ+uCgoKhDkMIIYQQQoxyGzZsaNBa5xz/+ohKngsKCli/fv1QhyGEEEIIIUY5pdS+3l6Xsg0hhBBCCCEMkuRZCCGEEEIIgyR5FkIIIYQQwiBJnoUQQgghhDBIkmchhBBCCCEMkuRZCCGEEEIIgyR5FkIIIYQQwqAR1ed5T30Hn330w5N6zQ59gBZTKS7ySdUlqJH1VyaEGKM0mna1lXb2YCGJdL0QKxlDHZYQQhjio5YWVU5Ya8Zz4VCH041kgjG41T6O2P7FEaWpD08kL/BlnHrKUIclhBB9CtDEQetjeMx7ul5r0M+TG7yKjNC5KNQQRieEEH3ThKmz/ING85ugwphD2eT6L0Cp4XPfUlrroY7BsCVLluih2GGwxeNhedlrvLDvIfxhD49d8BhFOUUnPQ4hhIjlQOsBvvzGl2n1t3N5/g18ueRqKmr38cTWB9naso6vl3ydmxfcPNRhCiFED6FwiB+t/hH/3PNPlmZ/gi/Pv5Flk6dhMg1NlbFSaoPWesnxr0vNswFpTiefnHkR35z9IKm2dL71zrdo9DQOdVhCCNGNN+jltvduwxv08q05v+TqWZ8mPz2Nc6bO52uzfsb5+Zfzu02/4+XdLw91qEII0cNjFY/xzz3/5NPTbuC6md9n6aShS5z7M/wiGqbyMpxkO3O5veQXtPnbuOuDuxhJo/ZCiNHvVxt+xc7mnXxnwd1McE5nUmYSAE6bmdxUJ1fk38LicYv5xdpfUNNRM8TRCiHEUZvqN/Hopke5uOBSzsj+LHnpDizm4ZmmDs+ohiGr2URWsg2HzuPWRbey8uBK3jvw3lCHJYQQAGxv2s6z25/ls7M/yyTnQlIcFpLtR5e1jE9zEA4pbl/4Y0I6xC/W/mIIoxVCiKPCOszP1vyMbGc2X5nzHbSG8WnOoQ6rT5I8D8C4FAf+YJhLp1zN1LSp/HLDLwmEA0MdlhBCcN9H95FmS+PrxTfT6gmQnWLvdjzTZcNkArsax1eLvso7B96hrK5siKIVQoijXt3zKlubtvLtxd+mw2clyW7u9uF/uJHkeQAyXTYAWr2a7yz+DlWtVby6+9UhjkoIMdZtqN3Aupp13FB8A6GgA60hq/N+FWU2KdKcNhrb/Vw791qyHFk8WPrgEEUshBARwXCQR8oeYX7WfC6a8gla3AGyXPbY3ziEJHkeAJvFhMtuodnt5+z8s5mdMZsnKp8grMNDHZoQYgz7ffnvyXRkcvWsq2lo92M2K1Id1h7nZblsdPiCmLDzlaKvsL52PeX15UMQsRBCRLxZ9SaH2g/xteKv0e4LEQprMpJ63r+GE0meByjTZaPFHUBr+HLhl9nbspeVB1YOdVhCiDFq95HdrD68mmvnXovT4qTZ7SczyYbJ1LMnalZyZDS6qcPPp2Z+imRrMn/e8ueTHbIQQgCgtWZ55XKmpk3lnEnn0OyOlMKmJ9lifOfQkuR5gDKSrITCmlZvgIsKLiI3KZe/7vjrUIclhBijntvxHBaThU/N/BTeQAiPP0RGHw+eZLsFi1nR4gngsrq4etbVrNi3QjpvCCGGRHlDOVubtnLt3GsxKRPNbj/JDgs2y/BOT4c0OqXUt5VSlUqpzUqpZ5RSjqGMx4jop6EWT6DrgfXBoQ841H5oiCMTQow1vpCPl3e/zPmTzyfLmUWrNzJqk+rsfaGNUooUh5VWT+S8z8z+DCEdkr7PQogh8fftfyfJksSl0y5Fa02LJ0D6MC/ZgCFMnpVSecAtwBKtdSFgBj43VPEYZbOYcFjNtHmDAFw540qUUjy/4/khjkwIMdas2LeCVn8rV8+6GoA2bxClIKWXeueoNKeFDn+QcFgzKWUSp44/lRd3vihrN4QQJ1Wbv403qt7gkmmX4LK6cPtDhEK61/Uaw81Qj4tbAKdSygIkAYeHOB5DUhyWrhGeCckTWDZhGf/a+y/ZNEUIcVI9v+P5rgQYoNUTwGW3YO6l3jkq1WElHIY2X+cAwMwrOdh+kA21G05KzEIIAfD2/rfxhrxcOeNKgK5ByRTH8G1RFzVkybPW+hBwP7AfqAZatNZvHn+eUuprSqn1Sqn19fX1JzvMXqU4LLh9IYKhyEjNxVMv5lD7ISoaKoY4MiHEWFHnrmND7QYun345JhW5lbd6gzFHbVKdkePR0o3zJ5+P0+Lktb2vndiAhRDiGK/vfZ285DyKsosAaPMGMJnAZZPkuU9KqQzgk8BUYCLgUkpde/x5WuvHtNZLtNZLcnJyTnaYvYpOibZ3jtycN/k8rCarPHyEECfNin0r0GgumnIRAB5/iEAwHHPUxmE1Y7WYumbPnBYnH8v/GG/vf5tgOHjC4xZCiCZvE2uq1/CJgk+gVGSmrNUbJNlu7bVT0HAzlGUbHwf2aq3rtdYB4AXg9CGMx7Dow6lrisGWwpl5Z/Jm1ZtSNyiEOCnerHqTGekzmJY+DYA2X+diQQP1gikOC+3eo4nyhVMupMnbJKUbQoiT4q19bxHSIS6eenHXa23ewIgo2YChTZ73A6cppZJU5GPH+cDWIYzHMIfVjO2YkRuIlG7UeeoorS0dwsiEEGNBnbuOjXUbubDgwq7XOnwhAJLs5pjfn2y34PaHutZpnJV/Fk6LkzerelTOCSFEwr1e9TpT06YyK2MWEJk5C4a0JM+xaK3XAs8BpUBFZyyPDVU8A+WyW7oeVgBn55+Nw+xgxb4VQxiVEGIsOL5kA6DDF4yUZJhj39ZddguhsMYTiNzDoqUbb+1/S0o3hBAnVJO3ifU167mo4KKuko0Of+S+k2yX5DkmrfVdWus5WutCrfV1WmvfUMYzEMl2S9d/bIAkaxKnTjiVlQdXStcNIcQJ9c7+d7qVbEBkDYbLwKgzQHLngpzoug2AC6ZcQJO3ibK6ssQGK4QQx1h1aBUazTmTzul6raPzXuSS5Hl0S7KZCYU03kD30edD7YfY27J3CCMTQoxmbf42SmtL+Vj+x7pe01rj9gcNj9pEk+xjZ89On3g6FmXh/UPvJzZgIYQ4xsoDK8lx5jA3c27Xax2+EDaLydDM2XAwMqIchqIPqY5jRm6iD7OVB1cOSUxCiNFvTfUagjrYLXn2BEKEw8ZHbSzmyGZPx96/UmwpLMpdJMmzEOKECYQCfHD4Az6W/7GuFpsQKdswOnM2HEjyPEhJvYzcjHeNZ3bGbEmehRAnzL8P/psUWwolOSVdr0U7ZwxkytNlN3cr2wA4K+8sdjbvpLq9OjHBCiHEMUrrSmkPtHf78A/RsrORUbIBkjwPmt0S6ZV6bN0zREafy+rKaPG1DFFkQojRKqzDrDq0ijMmnoHFdPRB0+GPfIh32YyP3LjsFtz+YLc1GtEHmow+CyFOhJUHV2Iz2Thtwmldr3kDkW25R8LmKFGSPMfBZes+7QmRh09Ih/jg8AdDFJUQYrTa1rSNBk8DZ+Wf1e11tz+IzWLCMoB6wSSbmXAYfMGjvemnpk0lLzlPkmchxAnx/sH3OWX8KSRZk7peG2mLBUGS57gk2SK9Uo9VlF1EijWFtdVrhygqIcRoterQKhSKMyae0e11jz9E0gBGnQGc1sj5x97DlFKclXcWa6vX4g/54w9YCCE61XTUUNVaxbKJy7q9Hr0HDfQeNpQkeY5Dks2MPxgmGDo6cmM2mTll/Cl8ePhDaVknhEiotdVrmZM5hyxnVrfX3f4QzgE+eKKjPO7jSs+WTVyGJ+ihvL48vmCFEOIYa6rXAHQr2YDIgmezSeGwSvI8JkQ/JXkC3UefT5t4Goc7DnOw7eBQhCWEGIV8IR9ldWWcMv6Ubq8HQ2H8wTBJA6wXtFtMmEyRUetjLRm/BJMysbZGZs+EEImzpnoNmY5MZmbM7Pa6xx8aUYkzSPIcF0c0eT7u4RP9VPVh9YcnPSYhxOi0qW4T/rCfpROWdns9+uF9oFOeSkVGeo4vPUu1pTIvc56UngkhEkZrzdrqtSydsLRbizqIzJyNpJINkOQ5Lkm91AwCFKQWkJuUKw8fIUTCrK1Zi1mZWTRuUbfXPXHUC7p6WbcBsHTCUirqK3AH3IMLVgghjrH7yG4aPA0sm9C93llrjScQlOR5LLGYTdgsph4PH6UUp004jbU1awnrcB/fLYQQxq2rXsf87Pkk25K7vR5tU+ccxLSn02bGEwj2WJ9x6oRTCeogG2o3DD5gIYTo1Fe9sy8YJhxmwGs2hpokz3GKPnyOt3TCUlp8LWxv2j4EUQkhRpOOQAebGzZz6vhTexzz+EMDblMX1Vu7OoCF4xZiNVll9kwIkRBra9YyOWUyE5IndHvdE8eH/6EkyXOcnFYzHn/P0eXoop7SutKTHZIQYpQprS0lqIO9Js9u/+CnPKOLDI+fPXNanCwYt4B1NesG9b5CCBEV1mFKa0tZMn5Jj2PurjUbI6fHM0jyHLckmxlvIEQ43H3ac7xrPHnJeTLtKYSI20c1H2E1WVkwbkGPY57AwNvURfXVMQhgce5itjVto83fNqj3FkIIgF1HdtHqb2Vx7uIexzz+IEqBwzqy0tGRFe0wFP201NvDZ9G4RWyo3SD9noUQcVlbs5binGKcFme314OhML7AwNvURUXb1bl9PUvPFo1bhEazqX7ToN5bCCEgMnMG9FjsDODxh3FazSilTnZYcZHkOU697dIVtTh3MU3eJqpaq05yVEKI0aLd3862pm09+jvD4NvURfXVrg6gJKcEszJ3PfiEEGIwNtRuIDcpl7zkvB7H3P7giFssCJI8x83ZR69noGuKQko3hBCDVV5fTliHex+16Uye43n4JNksvc6cJVmTmJs5V+5fQohB01qzoXYDi3IX9Tq67A6ERly9Mwxx8qyUSldKPaeU2qaU2qqUWhb7u4aXyCp3hfvYjhtag6+NKa48Mh2ZMnIjhBi00rpSTMpEcU5xj2PezsXKDsvgk2eH1YT3+OQ54IWAh0W5i9jcsBl/yD/o9xdCjF0H2w5S76lnSW7PxYL+YJhQSI+4ThsAQ53uPwi8rrW+WillA5KGOJ5BiXTc6Hz4hENwcD14mlAmK4uzi2TkRggxaGV1ZczOmI3L6upxzBsMYTYpbJbBj4M4rWaCIU0gFMZqNkHTXqjfDmgWOcbzp7CfzQ2bWZTbc+RbCCH6s752PdBXvXP8M2dDZchGnpVSqcDHgMcBtNZ+rfWRoYonHg6rGW+gs11d/XbwNEHWTDBbWWzN5HDHYarbq4c2SCHEiBMIByhvKGfhuIW9HvcGQtjjXKUeHfXxBkLgaYb6beDKgdSJLFKRBYrSclMIMRildaWk29OZlj6txzFvMJI8j7ROGzC0ZRvTgHrgCaXURqXUH5RSPYZWlFJfU0qtV0qtr6+vP/lRGuDsbFdHwANH9kP6ZMieARNKWJw8BYANdTL6LIQYmO1N2/EEPSzM7T159vhDcU952q3HtKtr2AlmG0wogfHFZCSPZ5pzvMyeCSEGZUPtBhaOW4hJ9Uw3o+VijhFYtjGUybMFWAT8n9Z6IdAB3HH8SVrrx7TWS7TWS3Jyck52jIY4LGZCYY2/cV/khczOT1jOdGaOKybF7KC0Rh4+QoiB2Vi3EYCFOX2MPAfDcT94osm3v6MN3I2QUQBmCygFWTNZlDKFstpSQuGeiwqFEKIv9e56DrQd6LW/M4A3EMZiVpFysRFmKCM+CBzUWkf3f32OSDI94jhskb9Gf/NBcGWD9WgvVnN6AUXJk9kk055CiAHaWLeRvOQ8cl25PY6FwppAApJnm8WE2aQIHDkIKEjLP3oweRyL0mfSHnSz88jOuK4jhBhbyhvKgUjby954AqEROeoMQ5g8a61rgANKqdmdL50PbBmqeOLhtJoxB9oJejsg+biHXFIWJalT2dWyl45Ax9AEKIQYcbTWlNaW9lvvDCRkpbrDaibUUg1JmWCxHz2gFIsmnAbAJukaJIQYgIr6CizKwpzMOb0e90ryPGjfAp5WSpUDC4CfD3E8g+KwmrH5mvCHwpGFNscymSgZt5Awms115UMToBBixDnYdpBGb2OfybMnkLjFNk7lJ+Rt73n/AibmFJJlTaa8+qO4ryOEGDsqGiqYlTkLh8XR6/FI8jzUaejgDGnUWuuyznrmYq31f2itm4cynsGymk04Ak14lROsPX9JCvPOAGBT9ZqTHZoQYoTaWB+pd14wbkGvxxO52CY53Io/FIKkrB7HVFImxSkFlDdujvs6QoixIRQOsblhM0XZRb0eD4TCBEdoj2cY+pHn0SEcxhVqxW3N7PVwWtpkpjlz2VRXdpIDE0KMVOX15bisLqanTe/1uDcQRimwx9HjOSop2EJQWwhYevaSRilKsgqpctdwxDMixzeEECfZnpY9uIPuXjd3gpHdaQMkeU4M7xFsZkWHObX34yYTJekzKW/ejtb65MYmhBiRyuvLKcwuxGzq/eESrRfsbcvbgXIEWwjY0/EGw70eLx4fWS1fXr221+NCCHGsioYKgD5HnqN7Y8SzO+pQkuQ5EbxHsJlNtKnkPk8pzi7mSLCD/c27TmJgQoiRyBP0sKN5B8XZvY/aQALrBQNe7NpLwJZ+dLOn48yfsBQTivIaqXsWQsRWXl9Oii2FKalTej3eNfJsG5lp6MiMerjxtmC1OwliJRDq/eFTMrFzxfrhD05mZEKIEWhr41ZCOtTnlCcksM2TtwWb2UTAltb1QDtekjODWa48NjVWxn89IcSoV9FQQVF2Ua+bo0AkeTaZwDYCezyDJM+J4W3F6sqI/M8+Hj7TxxWTbHZI3bMQIqby+khnnr6mPMNhjT8BPZ4B8LViNZvRtpSuDh69KcmcQ0XLLkKhYPzXFEKMWu6Am11HdvV5/4JI2YbDkpiys6EgyXO8QgEIuLEmR5Lnvh4+JpOZwtRplDdvP5nRCSFGoPKGcvKS88hy9ux+AeALhtE6QYttvK1gc2G3W/v88A9QPG4BHSEfexuk64YQom+VjZWEdbjfmTNvMIR9hC4WBEme4+dtAcAeHXn29162AVCSNY8d7Qdx+9pOSmhCiJGpvL6833pnTwI3SMHXAvYUnFYzHn8/yfP4UwHYJHXPQoh+RBcLFmYX9nmOxx8asW3qQJLn+HUmz7bkDMwmhTfYz7Rn7qLIZimyYl0I0Yc6dx217tr+R20StUFK0A9BHzhScdrM/ZZtTMmaQ5olqaukRAghelNRX0Fech6Zjt7b9x4tOxu5KejIjXy48LdHtrM1W7FbTf1Pe3Zuc1teu+FkRSeEGGEq6jtbPOX0Vy/YmTzH2+bJ1xr5056Gw2ImGNIE+1j0rJSiOG2GlJ4JIfpV3tD/zFl0kHGk9ngGSZ7j5+8AWwpAzGnPNFcOBc5cNknNoBCiD5saNmExWZiTOafPczyBEHarCZMpzsU20RIyewpOm7nrvftSnDWP3e4a2mSzFCFEL2o7aqlz1/X74T+aJ0nZxlilNfjawRbZlcthNfe5yUBUcfoMylt2y2YpQoheVdRXMDdzLnazvc9zvIEEddrwt4PZChZb1yh2X72eAUrGLUKjqahZF/+1hRCjzubOwcF+O2105kky8jxWBb2gQ2CPbI7itJoJBMN9TnsCFGcV0hRo43Dr/pMVpRBihAiGg1Q2Vvb74AHwBUKJ2ZnL3wG2yP3L3ll/2F/pWeGEU1EoKmo3xn9tIcSoU95QjsVkYW7W3D7Pid5j7JaRm4KO3MiHA1975M/Oh0/0U1R/o89FuYsAqJAV60KI4+w+shtP0NPvYkGtNd5gCGcidubytx9Nni0mTKb+k+eUpCymOsdRLpulCCF6UdFQweyM2TFmzhJUdjaEJHmOh7+zXtB2dOQZ+n/4zBxXgl1Z2CQjN0KI45Q3RDpZ9LfYxhcMEw6DPd6R56A/0qe+s+xMKYXDYu63bAOgKH0mFS27pPRMCNFNKByisiH2zJk3MLLb1IEkz/HxtYPZBhYbcHTas79Fg1abi/kpk6lo2npSQhRCjBzl9eVk2DPIT8nv8xxfIEH1gv7uM2cAdqu533abAMWZc2kOtHOw9UB81xdCjCq7W3bjDrr7nTmDBK7ZGEKSPMcj4O4atYGj056+GA+fovRZbG3dSyAUONERCiFGkIr6CopyivrdsrZrgxRbvMlzR+TPY+5hsToGARSPk9IzIURPXW02+xl51lrjDYRGdI9nkOQ5PgE3WJO6/q9SCruRac+s+fh1kO2NW050hEKIEaLd386elj2GpjwBHPEutvG3gzKB1dn1ksNqwh8MEw73XZIxI6cQp8lKRZ2UngkhjqpoqCDVlsqU1Cl9nuMLhtF6ZHfaALDEOkEplQ58ESg49nyt9S2JCEApZQbWA4e01pcl4j1PinAosjPXMaM2EPmF6K9PKkR2GgQor1lP4biSExaiEGLkqGysRKP7rXeGyAYDFrPCYo43ee6cOTtmlPvooucQSbbeHw8WRxpzXZMolw//QohjlDeUU5Td/8zZ0d1RR3bybOTu+y8iiXMFsOGYr0S5FRh5BcABd+TPY0ZtIDJy09+CQYDc9OnkWFOoqN90oqITQowwFQ2RKc/52fP7Pc/jT9Bim4C7x/0r+r79lm4oRXH6DLa27sUf8scfhxBixHMH3Ow+srvfzVHgaB/5kZ48xxx5Bhxa6++ciIsrpfKBS4GfASfkGieMP5o8J3V72WE14wtEpj37asOi7C6KUqZQLosGhRCdKuormJI6hTR7Wr/neQNhkuKtd9Y6kjy7crq9bKTdJkBR1jwCB95ke9P2mA9LIcToV9lYSViHDZedjYVuG39WSt2glJqglMqMfiXo+r8Gvg/0f6cejgK9J8/RXwhfrIdP+iz2u2s44j1yQsITQowsmxs2x3zwAJ09nuNtU+cDHe4x8my3mFCq/3abAMXjFgJQXpvISUghxEhVXh9psxnrHuYJhLBaTJhHcI9nMJY8+4H/BT7kaMnG+ngvrJS6DKjTWvd791VKfU0ptV4ptb6+vj7eyyZOwA0ma1ebuqjoyE3MuuesyNRsdKpWCDF21XbUUuepozC7sN/zAqEwoZCOf3fB6If/49ZsmEwKm8UUs+PG+PSpjLOmUlFXFl8cQohRoaKhgkkpk8hwZPR7njcQin+x8zBg5Cf4DjBDa12gtZ7a+TUtAdc+A7hCKVUFPAucp5R66viTtNaPaa2XaK2X5OTkHH946Ph71gsCXe1XYo3czB9XggklIzdCCDY3bAaMjdoA8bd56mPNBkRmz2K128SWQlHyZNlpUAgBdLbZNDJzNgp6PIOx5LkScCf6wlrrH2qt87XWBcDngHe01tcm+jonTMANtqQeL0dHhGIlz0mucUxPyqWiThYNCjHWlTeUYzFZmJ05u9/zulaqx93j2Q0osPQ2ABC73SZmC0Vp0zngrqHZ2xxfLEKIEa2mo4Y6T13MzVGgc3fBeO9fw4CR5DkElCmlHlVK/Sb6daIDG9a0hoCnR70zRKY97VZTzLINbCkUJ0+monmbbHMrxBi3uWEzszNmYzfb+z3P6+9cqZ6Isg2rA0w9HwGR5DkU875UnDkHkNIzIca66D0g1sizPxgmFE5A2dkwYCR5/geRbhgfcGJa1aG1fm9E9XgOeADda/IMxkduitNm0BpoZ1/rvsTHKIQYEULhEJWNlYYXC5o765LjEnCD1dXrIYfVhNaxFz3PzynGhJKWm0KMcRX1FVhNVuZ0fqDuizeYoLKzYaDfVnWdG5hcMKLKKU6GgCfyZ1/Js8VMmzf21ttFWfNhZ2TKtiCtIIEBCiFGiqrWKjoCHYZavnkDIeyJePAEPJCc0uuhrnZ1gVC/tYlJrnHMSBoviwaFGOPKG8qZkzkHm9nW73kJKzsbBvq9C2utQ0COUqr/v5GxJtAR+bOXmmcAp82ENxh72nNa5mxcZjvlUvcsxJgVnfKM1WkDErRBSjgEIX+fH/6dXclzjNkzewpFyZGdBsN65HUbFULELxgOsqVxi7GZs0SVnQ0DRjZJqQJWK6VeBjqiL2qtf3mighr2Ah5QJrA4ej1st5gJhyPTnv2N3Jgd6RS6Jsm0pxBjWEV9BcnWZApSC2Ke6w2GSXVa47tg0Bv5s4/7l9F2m1iTKE4u4Pm6dexr3cfUtKnxxSWEGHF2H9mNJ+gxNnOWqLKzYcDIT3AYeLXz3JRjvsaugAcsduhj//bow8dncORmx5FdeKMPNCHEmFLRUMH87PmYVP+341BYE4jxgdyQrrKznp02AMwmhdViitkxCKUoypwLyKJBIcaq8obI5ijF2cY6bSSk7GwYiDnyrLW+G0Ap5dJad8Q6f0wIentt8RTV1es5GCKNfkaJrE6KU6cSPPwuW5u2srBz1y4hxNjgDXrZ2byTLxV+Kea5nkRtaxtNnvsYeQZwWAx0DAKmZc4kyWSnvL6cK6ZfEV9cQogRp6K+gnR7OpNSJsU81+Pvfx3FSBLzI4BSaplSaguwtfP/lyilHjnhkQ1nAU+kzVMfog+3WLt0RUZu5gFHt7YUQowd25q2EdRBg5sLJGiletBLpMdzP/cwmzn2yDNgdqRRlDyJCrl/CTEmVTRUUJhdiOpjJv5Y3mA4/g//w4SRu/CvgYuARgCt9SbgYycyqGFNawj6+h15tphNWMyqqy1Lf7JT8phoz5DkWYgxyOjOgnD0w3hCyjYstl57PEc5rObYZWcA9tRI6VnzDik9E2KM6Qh0sPvIbkMlGwkrOxsmDA1haK0PHPdS7KxwtAp6ifR47nvUBgz2egawd26WIsmzEGNOeUM5uSaEuUQAACAASURBVEm55CTlxDzXFwxhMoE97h7PvW/wdCyHxUworPHH6PUcWbcxmaAOsbVpa3xxCSFGlMqGSjTacJtNSEDZ2TBh5C58QCl1OqCVUjal1HfpLOEYk7rqBfseeYbIL4iRac/oosFqdw317voEBCiEGCk2N2w2NOoM4PGHcVjMhqZH+xX09FuyAeCwHV230S+zleL0mYCUngkx1kQXCxqaOUtU2dkwYeSn+DpwM5AHHAQWAN84kUENa9GpSQMjz0YW3GBPpTh5MnD0F1EIMfod8R7hQNsBQ/2dIZLIxr25QLTsrI9OG1FdG6XEWrcBZKfkM8GeKR03hBhjKuormJI6hTR7WsxzE1Z2NkwYSZ5na62v0Vrnaq3Hde42OPdEBzZsGRx5dlhNhEKaQCj2Nt1z0qZjUWYq6uXhI8RYsbkxUu9cnBO7XhA6V6rHu7lA0Ac6HHPk2fBGKRApPXNNkpFnIcYQrTUVDRWGZ868gQSVnQ0TRn6K3xp8bWwIesFkBXP/Xf6O3eI2FkdSFrNdeTJyI8QYUtFQgUIxL2tezHOj9cfOeEeeg9Eez/3XPFvNJswGFz1jT6YoZRLVHdU0eBrii08IMSLUumup99QPIHlOUNnZMNFnBqiUWgacTmR77u8ccygVGB3j7oMR8MYs2YDuu3SlOGLsCGZPpciVz8sNZYTCIcymsfvXK8RYsblhM9PTp+OyumKem7A2dV0bpBi4h1nMsdttQteiZ4jUPZ83+bx4IhRCjACbOndGNjxzFkhA2dkw0t+d2AYkE0mwj91ZsBW4+sSHNkwFPTFLNuDoQ85Yu6cUipMn4Q662d2yO94IhRDDnNaaivoK4/XOCd8gxdg9zNCiZ1syc12TIqVnMnsmxJhQXl+O3WxndsZsQ+d7AqFR02kD+hl51lqvBFYqpZZrrfedxJiGt4AXnBkxT7NbzJhMxso2jh25qaivYFbGrHijFEIMY4faD9HsazbeaSOQoMU2BsvOILJRSosnEPs9lcLhzGBmcr6s2xBijCivL2de1jys5hgz60AwFB5VPZ7BQM2zJM7HCAUhHIi52CbKYTHYccOaxGRnLmnWZBm5EWIMGMjmKJDAxTYGy84gcv8KhjTBWIueoXPRYD4VDRWEwmN3GwAhxoJAKMCWxi2GNkeByM6CMHp6PIPBTVJEp67FNrGnPAEcNoMbpSiFcqRQmFLQVUckhBi9yhsiU54zMmYYOt8bCGNPSI9nr/EP/9FFz7E2SoGu5NkddLOnZU88EQohhrltTdvwh/0D6hQEkjwnhFJqklLqXaXUVqVUpVLq1qGKxbBAZ4/nAYw8GyrbALCnUOLKY/eR3XQEOgYZoBBiJNhUt4n5WfOxmmJPeULnYptEPHiCXrDYDZ0aXbcxkM2eAJk9E2KUi+5JUZJTYuj8rgXPttEzXhvzJ1FKTVVK/VIp9YJS6uXoVwKuHQRu11rPBU4DblZKxe7ZNJS6NkgxOPJsNeEPhgmHdeyT7SkUufLRaCobKuMIUggxnPlCPrY0baFknLEHD0QePnGP2oTDEPIPeOTZaMeNKY5sUqwu6fcsxCi3qX4TuUm55LpyDZ1/tOxs9Iw8x141Av8AHgdeAQzM3xmjta4Gqjv/d5tSaiuRXQy3JOoaCRdNns3GRm6iPVm9wRBJthh/1fbUrpGb8oZyTp1w6qDDjGrxtfDmvjfZ0rgFszIzP2s+500+z9BuQEKMZb6Qj7f3vU1ZfRmBcIDpadP5+JSPM941Pu733tK4hWA4yIKcBYbOD4c1vkA4/jZ1IV/kT4PJs91iwmQCn5FezxY7JouDotRpCRt5DoVDrKlewweHP6A90E5+cj7nTDqHmRkzE/L+Qoxm5fXlvH/ofWo6ash2ZnN2/tmU5JQkpM9yeX254ZINSODM2TBiJHn2aq1/cyKDUEoVAAuBtSfyOnELesFsA5Oxh1h0NzCP30jynEKaJYkC18SEjNz8bfvf+HXpr2nzt5FmTyOsw/x1+1/534/+l5sW3MS1c68dNc3KhUik1/e+zn0f3Ue9p54kSxIOi4PnvM/xqw2/4qvFX+VrRV+Lqxf7prrIugbDU56dyWv8G6QMrOxMKYXdYnDdBnT2q5/E7w+8jjvgJinGRiz92dG8gztX30llYyV2s50UWwoNngZ+s/E3XDjlQn502o/IdGQO+v2FGK0OtB7gfz78H9bVrMOkTGQ5smj2NvOHij+wbMIyfnLGT+IaBGjwNHCo/RCfn/N5w9/j8Y/N5PlBpdRdwJuAL/qi1ro0EQEopZKB54HbtNatvRz/GvA1gMmTJyfikoMX9Bl+8MCxI88GHj5mK1gcFKVO5cOGCrTWg0putdbc+9G9PL31aZZOWMp3Fn+HuZmR3dQrGyt5qOwh7vvoPjbUbuCes+7BMYCfR4jRLKzD3LPuHp7Z9gxF2UX8/Kyfc0ruKZhNZva17uPhjQ/zSNkj7GzeyT1n3YPNbBvUdTbVb2JSyiSynFmGzo+WTSRka24wXPMMkdINQx2DILJoMGkCYR2msrGSU8afMoggYW31Wm5991acFic/P/PnXFhwIXaznUZPI89uf5Y/VvyRsroyHvn4I8zONNZjVoix4IPDH/Dtd7+NWZn5wSk/4MqZV+KyuugIdPDizhf5zcbf8Pl/fp7HL3ycaenTBnWN6OCe0Q//EMmBUp3G1neMFEaGUIuAG4B7gAc6v+5PxMWVUlYiifPTWusXejtHa/2Y1nqJ1npJTk5OIi47eANYqQ5Hpz09/qDBb4jsNNjgaaC6o3pQIf5u0+94euvTXDv3Wh79+KPMy5qHUgqlFIXZhfzfuQ/zvQXf4p3973DLO7fgC/liv6kQo1xYh7lz9Z08s+0Zrpt3HU9e/CSnTTita4R5SuoU7jv7Pr635Hus2LeCuz+8G60NrGU4jtaasvqyAT94IAEjz10bpBi/hxneKAW61m0Ag5492960nW+98y0muCbw7KXPcvn0y7F3lsllObO4ecHNPH3h4yjg+jeuZ2vj1kFdR4jR5v2D73Pz2zeTl5LH81c8z7Xzru3avdRldXHtvGv5yyV/QWvNV978CnXuukFdZ1P9JiwmC3My5xg6P9rjeTR12gBjyfOVwDSt9dla63M7v+Lef1VFhlUfB7ZqrX8Z7/udFANYqQ6RaU+HxYzbyIIbiIzcOCIF+NHVrAPx/sH3eWTTI1wx/Qq+f8r3e04ttx5G7X2PLzomcfe0q/mw+kPufv9Hg0oChBhNHtr4EC/tfombSm7ie0u+12cXjC/O/yLfKPkGL+9+mWe2PTPg6xzuOEyDp8FwvTNERp6VSkCP56APlAksxkfMnVYzvoDxRc8ZVheTXBMHVffc7m/n1ndvJcWawqMXPNpzMZK/A/avZc6Rap6cfT0uzHxjxY1Utw9uoEGI0aKyoZLbV97OzPSZLP/EciYkT+j1vBkZM/j9hb+nI9DBd1d+l0DYwCZIxymvL2dOxhzDs9YJ+/A/zBi5G28C0k/Atc8ArgPOU0qVdX5dcgKukxjhMISMb5AS5bSZja1WB7CnMCtpPHazbcA7dR3xHuGuD+5ietoMbi76Ab7jS0Waq6B6E1iTYMICriz6Mt+YfDGv7HuDv2x6dEDXEmI0eXvf2/y+4vdcNfMqbiq5qd9yKX8wzNXTv8TpE87kVxt+xYHWAwO6VlldGcCAOm14/JFOG4np8Wz8wz8cfeAZKt2wJQOKorTplNeXD/hD+QMbHqC6o5q7TrsXh+m4XVz9HbB/DfjaIGc2edM+zsMLv4M36OHbK75OIDTwJECI0aDV38rtK28n3Z7OIx9/hFRbap/nhsOaHPsUvrf4R2ys28ift/x5QNcKhoNUNlYObLFgosrOhhkjyXMusE0p9UYiW9VprVdprZXWulhrvaDz61/xvu8J07XYZmAPnySbBfcApj2tJgtz02YMeNrzoY0P0eht4j/ybmfLIQ+rdjawo7Yt8gBrr4e6rZA8DiYthdQJkD6JG8/8CedmlfC/m/6P8kMfDuh6QowGNR013PnBnczLmsePlv6o3wR1b0MHq3bVs+lAKxeMuxmFmbs//MmAksRN9ZtIsiQxI93Y5igAbn8QRyJGbQa4ZgMgyWrpjMHAPcxkAnsyxa586j311LprDV+ntLaU53Y8x7njr8bfPol1e5pYu6cRtz8YGbQ4+BGgYfJpkDkNUnKZOec/+EnxN6ls3cNvP/jJgH4uIUYDrTV3rr6T2o5a7j/7frKd2X2e2+IJ8OGeRtZXNZMRPo2SzDN4pOwRDrYdNHy9nc078QQ9Ays7G4U9nsFY8nwXkdKNn3O05vmBExnUsDTAHs9RTquZUEjjN7Jo0OYCZaIodSpbm7YanlLZ1byLv+34O6fnXMbpk0tYPCWD/Ewn+xvdVO6vR1dvAnsKTFjYrVOIyergZ+f+mhx7Kj9efSe+aE2kEGNAWIf5r1X/RSAc4L6P3YfV3PeCll117eyuayc72c6iKRkszp/KJ/L+k7U1a1h1aLXha5bVlVGUXYTFZGStdoQnECIpIclzHCPPRmfPHOkURUvPDA4AaK25d939pFmzuHr6V1gwOZ25E1PxBsNs2NeM93BlZIOqvMVgT+72vRcU/ydX55/HE3v+wdqqt4z/YEKMAi/uepG397/NbYtv63c0+IjbT+n+ZrSGovw0iiel8blpt6C14oH1Dxq+3sa6jcAAZ84CIcwmNap6PIOB5FlrvbK3r5MR3LAyyJHnAT18lAJHGsVJefhCPnY07zB0jf9d9yB2s5OvFd/ErNwUMlw25oxPZca4ZNoPVlLb0gETSnptsZfiyuZ/Tvkhezw1PLLm5wP62YQYyV7Y+QIf1XzEHafewZTUKX2eV9Pipaqhg7wMJ8X56WS6bEzNdvHVkmvItI/nvnW/JKxjfzh2B9zsaN4xoAePPxgmGNJdI8BxGcTIs81iwmxWxjtuONOZ4xyH1WQ1XPf82t4VbGmq4JMF17Ns6gSyk+3kpTtZNDkdk7uBqt07CGdMA2dGr9///TN/yhRnDnev/Rlev+zOKsaGBk8D96+/nyW5S7hu3nV9nucPhqk41ILdbGJJQQa5qQ7GpTi4YPZszp3wKd7e/wbbGrcbuubGuo3kJuUy0TXRcJxuf2jU1TuDsR0G25RSrZ1fXqVUSCnVo6XcqDfAHqlRA6oZBHBmUOyITL0YqXve2ridD2re46L8T7NgYvdf6AJXgPE0UhXMplX3HfcZ0y/hU3nnsHzPS+zo7EF7IgTCAWo6aqjtqDWUbIixR2tNg6eBmo4avNF/cydAg6eBX274JUtyl3DljCv7PM8bCLGtppX0JCtzxqd0OzY5M5XPz7iBqradrNj7XsxrVjZWEtKhAU15Rj90x/3wCQVAhwZ8/wJIspoj5RNGONKxmSzMTZtmaORZa83vyh5lnCOfry38DGbT0bKZFLuF+bYaOrSVPbrvncyc9hR+fMoPOeBt4LGPTtykqNaaRk8j1e3VeIIySyd61xHo4HD7YVp8LSf0Ovesuwdf0Mddy+7CpPpO5bZWtxIIhSnKT+vWa9luMfPNRV/Fbk7ilx89HPN6WmtKa0tZNG7RgNZfuH3BxMycDTMxhzO01t2eGEqp/wDi3/5upAn6QJkj/ZgHINqeZSAPnwnWNLLsmZTXl/O5OZ/r9/SHNjyG3eTkpkX/2fMXumk3+Tlp7A9PY3tNG6cU9L2pwHeW/Zi3/7GOe9f8jD9c/teEbaASCod4vep1Xtz5IhvrNuIP+wFwWpycNuE0rpp5FR/L/5hs2DLGldWV8cy2Z1h1aBWt/shn8+iumJ+c8Uk+OeOTXS3LEuGB9Q/gCXr48bIf9/u7t6uuHa1h3sTUXs+7rviTPLPzMR4r/wMXTeu/CdGm+oFtjgLgDkTuG3E/fAb54T9ybQttXoML8mwuMFkpSpnK84feIxAO9Nm5BOCdfavY27aDmwt/SJrzuNhaD5Nu8pKSP599TV7Gp7tItvf+yFo69QIu33E6T+x+gcvmfJZpWYnr/7yjeQdPb32ad/e/S7OvGQCTMjE7YzaXTruUT838FCm2lBjvIkazBk8Dz257ljeq3qCqtarr9fGu8Vww5QK+MOcL5KfkJ+x666rX8UbVG3xzwTcpSCvoO652H/VtPmaMSybF0fPfYUFmDhdN+hQvVz3F7uZ9TM/oewbuYPtB6jx1LMpdZDhOrTWeQIhxqYm7dw8XA67g1lr/A4i7Vd2IM4h6QSBS62M1GW9X50xHKUVx+syY057VbfWsrnmH8/IvIy/1uIUCvnZor8OSWcD08em0uAPUtPQ9kpfmGsc3Z1/DuuatvL37FWOxxrCtcRtXv/IZ7nj/Dg60HeJzcz7HXcvu4r+X/jdXTL+CzQ2b+eY73+S6165jV/OuhFxTjCz17npue/c2rnvtOlYdWsW5k87ljlPv4O7T7+b6+dfjDnj5f2v+H1e8+Ek+TNCi1m1N23h1z6t8cd4XmZbW90YB0X8zkzKT+twh1GWzc9WMa9jRUsGaQ/3vG7WhdgPT06aTZk8zHGv0vhF3j9RBbJAS5bRFNkox1K6us/RsgSsPb8jLtsZt/Z7+eMVyUq1ZXFf0qZ4Hm/eCLZmCqTMwmxQ7a9v6fa/bT/svHCYrv1ybmPIzd8DNXavv5qqXr+K1va9z6vjT+MEpP+Du0+/mhqIbsJqt3L/+fi554RJe2vWStPwcg0LhEE9sfoKLn7+Y31f8nlxXLrcsvIWfnP4Tbl98O7PT5/LM1me44h9X8OCG3xIMGxxE64fWmgdLHyQ3KZcvFX6p3/N21LaRZDMzObPv3T5vKLkOkzLzaNkT/V43Wu88kOTZEwihNThj7bA8AsX8iZRSx97VTMASYOzdJQZRLxiVZDMb32jAYgdrEiUpk3m3di1N3qY+t6H9c+XfCekg/1n4hZ4Hm/dGerqmT2Gi2cqBJjd7GtrJTbX3OdJ2dckN/G3vq9y/4VecWXBhXLsPPrftFX627i6SLKlcN+1HLMg8G5vFzNQsF5Mzk1BK8YNTf8DLu17mwdIH+cyrn+Hu0+/m8umXD/qaYmRZfWg13/v39/AFfdy66Fa+MOcLXVs6Hz7iYTftzJ35GXa0lvLC/oe58a0buX7eTdy65OtxzVT8pvQ3pNhSuL7w+n7P293Qjs1ioiCr/22mv1h0NU9ue4S/bPk7p+X1/mAJhoNsrNvIZdMuG1Cs0W1tTaY4Z2YGsUFKVJLNjNaRbcL7+hDRjTOdxc5ICdmG2g0U5RT1etqupv1sbvqIz874Ci7bcXF1NETa0o0vwmYxMTXbxc7ado64/aQn9d6nOittCjdMv5Jf7XiGtQfeZ+mkswb0cx5rf8thblzxdQ51VHF27lVcMPEaXJZUxqXamZWbgsNq5psLv0llQyX3fXQf/736v1l9eDV3Lbura2MKMbo1eBq4/b3bKa0r5ZxJ53D74tu7RoE7fEG21bSSb7qYc7LrefXAH/jD5sdYfWgtv7vgITKdg+/++86BdyhvKOfu0+/udzaurs2H2xeiOD+t3/vHlPQJnDH+At4++E/afN8jxd77729pbSkptpQBdQqKlp25xmLZBnBsNhMEqoBPnpBohrOgFxyD+4V3Wi00tA9gJ7/jHj4XTLmgZzihIP+seoG56QuZn3PcL3PAC62HIS0fLDYUMC3bRfnBFmpbfYxP6/0BarE6uaPkG3zlwx/z1KbH+OriW4zHfIzfl/2F3266hxmpxfz09PuYkj4ObyDEgWYPO2vbaerwU5yfjtVk5apZV3HuhNP53qof8l+r/os9LXu4ZeEtJ6SMY1fTPl7c+Qpl9Rs43HGAoPaT6chgdsZszsw/k49P/nhX8nay+EN+2n1uXDYn9kGMCsYjGA7y3oH3ePfAu2xr2ka9uwGTMjPOOZGirAVcPuMSFuTOOyHX/vuOv/OzNT9jRvoM7j/7/q6HjtaardVtHD7iIT3JyuzxKZxuu5DLZi/jJ2vu5vEtj1DvbuanH/vhoH5HSmtLef/Q+9y66NZ+R4BbPAGa2v3MzE3GYu5/gi7dkcKZEy5gVfUKmjx3kOns+b7bm7fTEehgce7iAcXrCSRosU3XyPPgkmeIPAgNJc+OdHJsKUxJzmd97fo+R8ee3vI3QHHN/Kt7HmzaGxlISIncB/MzkqhqdLOnoYNFk/ve5OWaBTfx16rXuf+j+/hr/hn91oL2pbJ+Jze9dSPuYAd3nfJrLpgWScLrWn0caHKz1t3Egvx00pKszM+ezx8//ih/3Pw4D1f8nv2t+3nk44/0OeARj3ZfBy/tfINVh99nb+sOOoJtOC1OpqROZknuEi6ZegmTUicl/Lr90VrT6mtHKUWKzXXSy+92H9nNq3tepbS2lIPtBwmEgqRY05meNptzJp3LZTPOx2Y2vimQUfta93Hjihtp8jbx0zN+yhXTr+j62evbfFQcOoJJKWbmJpOVnMUFs/6XZ7f+g/+r/AWff/U6/nTxE+Qm991Wri+hcIjflv6WgtQCrph+Rb/n7m3oIMluJicl9nPls3Ou5t/Vr/H3ba9yfclnez2ntK6UheMWDujflDtRazaGISM1z18+GYEMe3GMPDttZvzBMKGw7rYgpk+OdOY7x+M0O/io5qNek+fX9qykyVfLN0tu6/n9R/aB1pBR0PVSToodl90Sc/T51BmXc/a2Z/jjtqf5dOF/DmiKGeCpzS/wm02/oCRrGb+74Nck2yPJqMtuISvZzsFmN9tr2ig70MzC1HZMR6rI9Hfwu4Kr+TlW/lDxB/wBN9899Y6E3Yj3tRzkvrW/4f3qN9CEyU+aQYGrCJvJgTd8hHU1H/Fa1Wvca7uXG4pu4PNzP5/Q+trjbWncwjNbn2NN9Rpq3YfQRBZPZtlzWZy7hE/P+Q+Wjl96wh5EYR3m+Z3P84fyP3C44zCptlSmpsxjTupMQjpInfcAz+36M3/b9SRFmafy42U/ZG628dGGWJ7e+jT3rLuHM/PO5P6z7+82Uld5uJWaFi8F2S6m5xx9GLvsGTx60QP8cOVPebnqGVLtafzgtJsHdN3odGe2M5tr5l7T77n7GjuwmBV56cZaU35h7md459Ar/HXLy9y0uOfK9/U16wEGnDy7/SHGGXj4xRT0RtZr9NJxJxZH17qNEFmGviFyz1iSMZc3a9YQCod67HbqC/p5c/+rlGQtpSA9r/v3e1vA3QDZs7riNZsUUzKT2FXXTosnQJqz9zpquzODW2d/gR9UPMwrO1/kk7OuGtDPeqClmpve+jpBHeR35z/OkomFXcdSHVYmpDkoO3CE0gPNLMq1kNa+F3NHPTckTWPWrOu4fefTfOlf1/HHi5/st+fuQATCAR7e8ATP7niSjmAradYsJrvmMiUpnaD2Ud+xn4fLHubhsoe5ZNolfKPkG0xOnZyQa/em2dvM37Y/z1tV77KndQf+cKQc0GF2MSdjPpdOv4grpl92QgcitjVt4/7197O2ei0WZWF62hymJi3EpKy0+OtZW7OKdw+9xoMbc7mh8Caumf+phN1P9xzZw5ff+DJaax6/8PFuMyv1bT7KDx4hxWGl+LgFejcu/jRTMyZyx6pb+OobN/H0pU+Q6kju7RJ9enXPq+xu2c0DZz/Qb7vLxnYf7d5gn2s1jnfWpFPIdU7ilT0v9Zo8N3mb2NuyN2bCfjy3P4TZPPra1IGxso37gJ8CHuB1oAS4TWv91AmObfgI+kGHB1UvCEdHbtz+YK9F+z04M7CazCzInMf62vW9nvLizhdJsWZwxcwLux8IBeHIgciGKLajSYlSimk5LioOtlDX5iM3tY8PAiYz3yq+kU+vvI0nyv6P25beYehnBHj/wEfcX/pT5qYv5PGLHsJu7fmJPz8jCQth9m/+gAPmVqbkTYScOVjNNu5Mn4Kt3MKftv0Fu1bccprxa/dGa83ftr3I/RvuJaRDXFbwGW4o+RJT0/PwB8PUtnqpauzAHwzht+zmpX1/5oEND/Dczue496x7mZ89P67rH29L4xbuWXcvG+tKsZkczEpbyBnjzyfTkU6rv52dzTtZdejfvLn/n8xMn8V3l9zO6XmnJzSGvS17uXP1nZTVl1GSU8KNRd8mTZdgUmYKslxMSHPgsJqpbW/iz5V/4287n+Tz//o01839Ct9e8o1BjeQd64WdL3DPuns4b9J53H/O/d0Wk+1vdFPT4mX6uGSmZvecOrSYzdx37p20vdnKU9t/x+SUPD4//z8MX/v9Q+9TWlfKj5b+CKel76TY7Q9S1+qjINsVc9Q56tSJJeS7pvPK3hd6TZ431G5gcspkxiWNMxxvIBQmEAwnqMfz4D/8O6xmzKYBtKuz2MCewuLkyTwfWMHOIzuZkzmn2yn/2v0urYFGrprZS3LbtBdMFkjvngDmZzipauygqqGDkkl9zwJePO8a/rznZR7a+BAXT7/c8Mhju6+Dr6+4mY5gK49+/I8smVDY4xyX3cLiKRlsrtzMgY3bcU5Ix5Y1FewpnJ09i9/Z0/jG5kf4+mtf4olL/0KKve/d3ozY3rST7753B1VtOyjOXMpXCr/C2ZMjH6yb3X6qGt00d/gxWVvY2PIqf93+LCuqVnDb4tu4Zu41cf97PVabv42Hyx7m79ufwx/2MSlpFudOvJyJKeMJ6zAHWg+xuamUn639f/x6wy+5vuh6vjT/Swkd+fWH/DxY+iBPbX2KdHs6N5fcwkzn+eiQi3GpdqZkukh1WvCHgry2513+uPkP3Lvhf3it6p88cO7PGe8aH9f1D7cf5oYVN6BQLL94OVPTpnYda/MG2HyohRSHlUWT03u9d1w47Qy8wXv48Yff5RtvfZsnL/ldjw+W/f3sD5c9zLyseb0OqB2rqrEDu9XE+L6e88dRSnFpwSf549aH2Fy/k8Kcmd2OR+udB/7hPxj/eo3hSmvd7xdQ1vnnlcCTQCawKdb3nYivxYsX6yHhadF627+0bq0e1Le3evx6RWWNrmnx1xcS/wAAIABJREFUGPuGcFjrnSv0o6vu1oXLC3Wzp7nb4Ub3Eb3gyYX6++/c3fN7G/dEYnU39zgUDof16p31es3uhv6vHwrq77/6n3rJnxbpuo46QyHXtNXqZU+fqc979hO6rr2p7xODfq2rPtDV617Q76/9SO+tb+8eo7dN3/XG13Xh8kL9QvkThq7dm0AooH+48se6cHmhvvKFa/T2hqpez/MHQ3rj/ma9orJG76lv16sPrtbn/+18veBPC/TTW54e9PW7X8Ov71t3ny5+slgve/pM/d/vPKzLDh3WoVC423nhcFjvqmvSP1u5XJ/77EW6cHmhvuPfd+hWX2tC4nhn3zt66dNL9RnPnKFf2vWSPtjUod/aUqPX7mnUHn+w1+852FKnr//nbbpweaH+6utfjyuW1QdX6+Ini/WNK27UvqCv27Ejbr9+e2uNLtvf8/f2eB6/V1/14rV64ZOL9Oa6bYauHQqH9FUvXaU/8dwntD/o7/fcLYdb9Ntba7Q30PvfSV9++9ETunB5oS6v3drj2mc8c4a+c/WdA3q/I+7IfaO21eB9oz97V2l94KNBf/uHuxv0RgP/bbrUbtGHy5/RhcsL9VNbnupx+Mv/vFkve/os7Qsc99/C79Z622ta127p9W131bXpFZU1us0b6PfyH5T/SRcuL9RPV/7JcMg3v/ldXbS8WL+8463+T2zYpd0Vr+j1q97Q63ZWd/93HArpVRVP6wVPFusvvXS19gd8fb9PDO/se1cv/tMp+rSnztR/qXilz/P2Nx79d3yotUbf/NbNunB5ob5pxU26zdc26Osf67397+lz/3quLl5erL/y6nf1S5WlusPX87/BEbdPP7Vxpf78P27UhcsL9aUvXKo3129OSAzV7dX6C69+QRcuL9R3f3C3rm1v0qt31ut3ttXq2j6erYFgUD+47gm9+E+n6DP/crbeVLdp0Ndv8bXoy164TC/7yzK9rbH7fccXCOn3d9Trf++o6/Neeqz/2/BnXbi8UN+75reGr//Ulqd04fJCvfrQ6n7PO9IRuW/sa+gw/N5aa32otVYXP1mi/3vlL3ocu2/dfXrRnxb1uG/HsnpnvS4/cGRA3zPcAOt1L/mokY+l0aGhS4BntNZNJyCHH94GuUFKVLROsMNncKWtUpCUxRLnBCAyanWsV3a9SVAHuGLmpd2/LxyG5ipwZkIvCxKUUkzJdtHmDdLU4e/7+iYz3yy+gaAO8mjpb/4/e+cdJ2dV9v3vmd5ndnZne8umbXpIoRO6IA/iY1fEB0EFpAVUUINKU1FEeRArVsT31cf+6qMiSFNKAklIQirp2d7r9HK/f5ydLdnd2Zmd2exuON/PZz+T7Nxzz9l7d8657uv8rt814XATWoI7/rWBUDzI187+Fj772M0MSCSgcRuEeiiuPRN36XwOtvXTExiywRJmB184/2HOyFvEfa8/zGtHn53w/Y8nGAty0z9v4S+H/8jby6/iF2//KQvyx7bgMep1rCh3U+y2cLC1n2LTcn5/xe85u/RsHnj1AR567aGsPKmb/c1c849r+MXuX3B20eV8YcXj3H7qx1lRWjKqiEMIwVxfHree9mG+sPInXFp2FX8//CRX/vVKDnUfmvQYAH6x6xfc+tytVLmq+N07fseZRZewt7mPPLuJ1VV5I7YXh1Pm8vH9tz3Eh+fdxqstr/DRv19Ld6g74/c/3HOYz7zwGeZ55vGtc781IhsVjSfY2dCDSa9ncenEmTqL0czD5z+ExWDj0y/ckZYf9JOHn2Rf1z5uOuWmlJ0EQ9E4TT1BStzWjLca37XgMgQ6/rT/byO+v79rPz3hnsz1zrly2oBJuwUlsZky8HoGsOVTYnJTZiselKwk6Qn1sa39Fc4uuRCT4bjfRdcR+ThMcjacijwbep3gSHvqZiinz7uCNc4afrTjR2l5Mv/Pnj/yQuOTvG/uNbxj/oXjH9hdB+1vYvWWU7HsHHoigv2t/UPP63SctfRK7l+xns1de/navz474XuPxW/f/C3rn1uPz1LBjy/6FR9aOn6haYXXxrJyN/3hKPXtRr517iNsOG0DLze+zNVPXk2zv3lSYwA5t3/n9e9w87M3Y9O7Wb/4Ue5c8yXesWjlmPp3t9XEB5efw4Y1X+f6BQ8QiIT5yN8/wm/2/WbSYwDZRffKv17Jge4DfPPcb3Lnmrt4szFOJJ5gVUUeheNkWA16Pbeu/SjfPPvH6DDy0Sev4cWGFzN+/3gizp3/upP6/noeveBRFnpHWiHubuolEo+zosIz7lw6nOtWXskZhRfzy70/4sW6jRMe74/6eWzHY5xWfBpnlJyR8tgjA5KzUk9mO02lzkKWeU/jhYanR615m1s2s7RgaUa7CImElruajRlIOsHzX4QQe5EuG88IIXzA1HUvmIlk4ZEKUq9nMerTt6sDsPtYai3GojePkm48deRJ8s0lnHl8ZX9fkxyrdw7jUeKyYDLoONKRevGpKDuddxedzu8P/YW6vrqUxz6+85e83raJq+bfzNqyReMf2L5PahmLl4KziNoSJ2aDnp2NPcTiQx9Wg9HCQxf/gAqLj9tf/Dz1GdjYReIRbn32Vl5pepkP1tzOfes+g8OS+gMvhGBJqYtip4GGI/sJH9vPf9dew4eq3s7jux/n7pfvnlQAfaj7EFf97Sr2d+3nutq7eW/Ves6aU4HXnno8DrOBM+YU867qj3Hzom/QE+7lyr9dOSoISQdN0/j+9u/zjc3f4OKqi3n80sex6fN5o74Hh9nA8jL3hDp8k0HPp0+/hpsWfZVDPQe55h/X0hXqSnsMvZFebn32Vox6I49e8OgoLeTepj5C0ThLy1wY05RJVLiLuHPVvTT4D3P/yw+mPDaaiPKdbd9hft58LptzWcpj67sCaBpUTeCwMRZlrkJqPSt5oeGfI2zLkje/mQbP/khM3kdna/OUSEA8AimkKhNhMxkIRtK0qwN5A49gdd5CtrRsGXE9/nrgGaJahMvnXjryNfEo9NSDsxiMY4/VZNBRlmelpTeU0sFIWN3cPP8DtIe7+J89v0o51LreOr6x5QHmuZZz52kpdPShHmjdDbYCKF5OoctKZb6Nus7AqILwy1dcyzVz3slv6v7J77f/KOX7H88f9v+B+165j1r3Gh46+wcsKZrYH7jQaWFZmYdQTwsHdm3hA45avrfqThr76rnmyWsmFUBH41HueOEOfrjjh5xX+h/ctPDbnD9nNQuKnCl1tHqdnE/PrzqH9Yu+y/L8Ndy/8X6+tflbk7Lz29m+k4/+46MAPHHZE1xQcRHb6rqJJhKsqsrDbZtYCnnunOV86+yfUGip5NZn12ccQD+y9RFeaniJDadtGPU5rusM0N4XZp7PiSsdWSag0+n42rn34rOUctdLd9Ef6U95/BO7n6Az1Mn6VetTXvv+cIy2vjAVXlvakrPhXFR5MV2RVjY1bhv8Xk+4hz0dezi95PSMzhWKSZu6k7FBCqTXnvtzwBnAGk3TokCAt5rbRrJSPYsiMptZn37mGcBWgElnYEXewhHBc1ugg52dm1lXetHoD9GALyp237in1ekElV4bnf0RelM1PtDpuX7ZJzAIHd/f8t/jHlbXW8ej2x5hied0PrnqqvHP198ms0qeKukCgsz6Li1zDXRxG+nh6rJ6+e4F35FZ7eduJxqd2K0klohxxwt3sLFpIx+s+RQ3rb4qrSwAgOg+xuLgVkrDB2ioO0J/TyefLz6fG8ou4k8H/sRXXr43o4l/d8durn7yamKJGJ9d8W2WuM9hebk7rYkeZJHpKZV5LHCv4NNLv4fPWsgN/7yBf9f/O+0xAPz4jR/zvW3f44q5V/DgugfREka2HevGZNCxchxd3lgY9Tr+a+Ul3LjoqxztPcqN/7wprYyepmnc8/I91PfV8/B5D1PqGNkFs74rQEtviLk+x7gWZOPxzoXnc1HZe/jz4d/yauOWcY/74/4/UtdXx62n3JpSAxqNJ6jrClLotEw6YL2o8m20BOvY2rxr8HubWzZTYi+hzFGW4pWjCYTjWAf0xlkRn7zHcxK7WdrVpa171hvAmscaeyVd4S4O9QztnDx19B94TAWcVb5m5Gu6j0EiBt7xvbeBQc/aox2BlMetnncZZ7kX8JM3fow/OnayQNM0vvji/aDBvWd8BbNxnM9nPCZ3zfRGKFkxWMg4z+fAbjawu7GXSGzkDfb6s+7hTO8SvrL9u+yofynlWJM8efhJ7nn5Hmrda/nsqgdYXJxm0WGoB1/HZpYk3iTcdohjTc2caS7ksdpr6Q6287G/X02LvyW9cyF372557haeOvoUH629mStKb2N+oWfMWoSxEEKwuMTF3Pwirqy6h/+oeg8/2/Uz7nnlnowSEW92vcn1T1+Pw+jg8bc/zjz3fHY09OAPx1he5k47WAVYVV7BPWu/TaGlkvXP3jbYtGginjryFD/b9TM+sPADvG/B+0Y81x+Osb+1j3yHicoMb7i9NiefW3MPXeE2Htg4fmfMrlAXP9/1cy6svHBc28ckRzv86HWCirzJFWu+Y95F6IWBvx58cvB7m5s3o6FxWslpGZ0rmSx8ywbPAJqmdWmaFh/4t1/TtMnvA81GYiHQmyZVqZ7EYTYQiMTTD8CMloGim2r2de4bbPX5l/3/IEGCy+cel0Hrb5O+qN45UvaRgrI8K3q94Gh76sWnsPgUPlh8Nn89+vSYsgFN0/jSS/ehw8AdqzeMb4QeC0PzdjA7wTeycMhjM1FdYKe5JzSqiUtFQS33rrmTnX1HePjle1KOFeBbW77Fs3XP8q7Km7h62fvTC8YScWjYAq270Vk9VK68kHD1+WzVL6e/8jxuPOUWri27gN8c+APfePm+tH5/B7oOcP3T12Mz2Niw8ru4dNUsLnGT78gseJEBtAe30cfNi77FHHcN659bz6amTWm9/jf7fsO3X/82l9dczv1n3U8ioeP1ui4SmsbKCk/GsgSLUc/7l17If829i10du7jjhTuIJ1IHU7/b/zuePvo0t666dZS5fl8oypstfXgdpklleoUQfO6028kzFfKll+8hEh8tRQrFQvxw+w9Z4VvBueXnpjxfQ1eQeFyjqmDyLgHvnH8JOnT85YCUbmiaxpaWLRlnnUEuzLZxOuplRBYez0ns5gylZwCOQlbbpPQsuWvSHephW8cmziq+cGShVFJyZisAS2rpjsWop9htobE7OCpgHYHVw80LPkh3tI9f7vrFmIf878G/s6XtFd5b83GWF1ePf67W3fI6lqyUBZED6HSCpWUuYokEe5p6R7xErzfw4AWPUmhyc+eLG+gLplY8bm/bzl0v3kWNcym3LvkyK8t96TlE9DTAsY0QC5Ffswrnsss4YF/F0bzTWVb7Hn6w/FY6Qh18/O//RXdw4h2jaDzKp57/FC83vMz6FRtY4fxPSvOszCvMrItickev0GXjQt8NfHjhtfxh/x94YNMDac2jx3qPcf3T12MxWPjJJT+hwlnBvpY+Ovsj1Ja4Mp5PhRCcWlnOnSu/icuYz83P3EJ9X33K1zT2N3LPy/ewrGAZn107UoKTSGjsbOjBoNOlJTcbi4tqTuXC0vfw58O/Y1Pja2Me8+M3fixvZk65JeW5QtE4zT0hyvKsmAyTi1XybR6Wetfy78ZnB39HG5s2YjVYWVaQOnA/nuRckfXO2Qwld6W4JzNZVKonsZn0xBMa4VST/fE4S1hrK0VD9pQHeK7uOXyWMlYXH+e/23lQjtFZOsaJRmLU66jIs9LaF0qtY9TpuWbptVh0Rr43Rvb5fw/+jc2tG3n3nI+zqqx67HNoGjTtkIvjsIzNcGoK7HhsRvY09w7qPJNcXPs+PlT5Np448r88d2D8zod/2P8Hntj9BOuK3s2Vi65Mz2IsHoO6V+WNh68WKtZidHhZWeHBoBNsa+gn7KritnO+yofLzueJA7/jl69/L+Up63rruO7p6zDqjHx25cPo4gUsLHaO6609LpoGvU0427ezOvo6xY1buNP7fqpsJdz67K280Za6++RTR57iyxu/zLryddx31n1ommBHfTehaJyVFZ7BYChT3FYjH1jydt5ddRMv1L/At18fXxN/sPsgD776IGeUnMHVS64e8Vw8ofHGwMKzJE07pbEocrq5cdlnafAf4fuv/2TU87/a+ytag60TbnfGExrHOgN4HaaMslmjxuMoYHHeGv7V+AyaprG/ez+doU5OLT41o/NomkYwGsNhzqXHcxaZ52TdRibSM2cxFeZ8ii0FbGqWN3x/P/QscS3K22uOk2z01ElpSf7ctE5dlW8nntCo60qdAFg652LOz1vC47t+PpiASNIb6eXrr32dCvsCblp99ThnQHrm9zbIsdlG+zc7LUbm+Zy09YWpP248bruPr5/1ZZrD3dz3wp3jBo1N/U2sf3Y9LlMB19fex9rqovR2HDoOQvMOsObBnHXgqWBuoZMil4X9rX5a4g5WLP0w312zgcZAK+uf+jjhFDUC8UScu168ixcbXuT2Uz7PHPMF5DtMLCqeRHAY7ke07GJZeBvVXS/xtmAlHyq7lF/v+3XKeQOgI9jBdU9fRywR40cX/4gyRxmH2/00dAWpLrCnbSF5PDqd4Kw5Vdy0+AEi8Rg3/vMmeiO9Yx4bS8T43L8/R4IEX1/39VG1Egfa+gft4CZrxSaE4M7TbsNrKub+jV8Z1YGw2d/Mr/f+mnfUvIO5ntSfjeROTKpugulwYeXFdIZbeHVAurGpeROrilalrBUZi/5wDJNBN+lAfqZzcv5UuSYWkpngLLBnWjQI4CxmuaMSi87EpuZN+CN+dnZuYW3h2eiHb7f7OyDYJbc708yOl+fZEAKOdaZefLxFy7iq9Dz+Uf8c+zr3DX4/EA3w0OaHqLAt4IZTPjJ+YNJ1ROqcC2tl5nkMhBAsLXMjgDcaekbpKj9z1v0sclTwhU1fpqn32KjXb23Zyv0b76fWvZqrF97EgnQyJIkENL4udYylp4zQiVuMelZWeIglNLbVdRM3WLnzvIe4yLeKb7zxQ57d/6exf9RQF9c9fR2RRIQvrH6YRCSfOT47FZlOZhE/HHsFmrZBuA9nQSnl1Qsh7mBDwXvx6q188unrOdxzeMyXv9H2Bhte3MAK3woeOvchDMLArsYeugNRlpa6M5ZHHE+Jy8KVC9/PGb7L+enOn/LUkadGHROOh7njX3dgM9r46jlfHSWX2NvcSyAcZ2mZO2sP0PcuuohTvOv4xZ6fjtia7ov08ZOdP+Gs0rNYW7w25TmSWcw5+dl3h1tXdh5toUb2dRzklUbZUvyM0tRFPscj22HnKGuTZc0GDNVtZDR/Ga0Iq4fTPQvY1LSJeCLOc8eex2XM56zyU4aOSySk5MziGTM4HQuH2YDPaaauMzCiXmIUA9nn/liQx3eObD/80KuP0BvpZv2KDTgt49xYRPzQsksGp/nje51XeK14HSb2t/TTf9w1WlFxDjfXXsmTLZv4046fjnptKBbilmdvIRAN84kF93PmnKr05GZdR6D9TXCWQNkaKSlhKOPrsRnZ1dhDVzDG6tp385VTbmNr95t88ZlbxpROaJrG1179Gn8/8nduWHYrVaYLcFmNLC/3ZNbhMpGA1j1w5EXobURvdVI9txaDo4DLdKfyroJT+fEbP+aJnY+P+fJIPMLtz99OR7CD71/0fWo8NTT1BGVBt9vCXF92n1GzQc8l8xbzsfn3cKzvGBv+vWHMm5rHdjzG662v84XTv0CFc2Tzmfb+MMc6AlR4bRRkmAE/nhKXm6trb+Fo30F+tWdkYeX3t38fDY0bV96Y8hzhWJzG7iDFA1aj2XD53IvQCT1/PfQULf4WDvcc5vTizPTOAP5wfNJJmtnAhJGWEOKZdL53UhMNZZ95Hsgg+cMZZG5MdkxWL6vcc3ml8RWeOfpvYlqUCyrPH3lcxwGZVXKn313KYtRT7LLS2B0kHEsxJp2eq5d+FKfeyneHZZ9/tOPndIbb+cTS2/Dax7k2oR45uTsKR3m2jjWeRSUueoNRDrWPLJ4wmWw8dM7XiSXibHj+0yOkAp2hTj7zwmfwmor5xMJ7WFlRMPFEr2kyWzOsePF4nBYjy8vc+MMxdjT0gM7IVy98lCXOSj638X52NYyUTkTiEW577jZaA63cteohEuFiyvKszPVlZoJPfyscfQUiASheDjXnQfEy8mtWUrx0Hf2ec/lS9Y3oNY2bn7pu1DZsi7+F9c+tp8BawCMXPIJFb2F3Uy+tvWEWFDnHrUqfEE2TmtQjL8Gb/2B+z8vc5jmNGmsNX3jxCxw4rqjzm5u/yf6u/dx/1v2jmkU0dAdp6paNUCYqnkwHg17H7atvJ67F+drGhwe///OBbOOtq1J3ykwkNI50+PHYjOTlYDyXzJGfzycPP8srTa9Q7arO2F82GYDlZPGJhUHoRsgNJoPdrB8VGE6Is4QzHFX0RnrZ3rad19s3scp3Bgb9sAW+r1FKItLMOiepzrcTi2s0dqeuX19QfT6X5i/nl3v+Dx3BDgAOdB3k/x38HeuK38FFc8duqU4iAU3bASF3zVLsXCQDVr1O7vAcH9Bfs+bTnJa3iAd2fI/DbTtHPPfgaw+yr2sfV9V8ngtqlo/bAGYEvY0yQHUUjrmjp9OJQeeHbfXd9IWiXLrsam6r/Qh/b97Ioy/dN+qU/3fv/+XX+37NlQs/whL7O7EY9awo92SmuY/HoG7TQH1LhZy/ylZjLFnCgpVn4y8/h8t8H+N8zxK+seWbPH945I23pmnc+8q9vN76Ol8++8ssLVhKR3+YPU295NlNLC6Z/C4VwS5o2Ar7n8ZZ9xzvt0R5X/4VvFD/Ao/vGhnIb2nZwg93/JB31LyDy2tGOp2EonF2N/ZiNxuYX5jh/D4OVy69jPmulXxv+/cGd0gO9RziTwf+xAcWfmBUrcjxHO0IkNC0tDXpqfDZ86h1r2Rj84u82vwqQMZ6Z5AFz463YvAshLAIIbxAgRAiTwjhHfiqBibWBpwsJOKQiGYdPJsNeowGXeaLj6eSM5w1HOo5xP/s+w02vZN1lcMKbfrbINgJeXMy1mRXF9hIJGS1cCpcvsVcXXY+zzW+yM72nbT6W3li989Z6T2HyxecNfaL4tGhApui9LRSRS4LpR4rR9plEdlwKguXsWHZ9Wzu2svPtjwCSBulDf++i65QNx+d/0XOmFOW3hZR6x7pTFKwYLB4cSzyHWZqS1x09kfYVt+Nyejk0Yt/iMdo55YXPkVrj8yCa5rG/RvvZ2vrVm5ZfhfG2BxKPBZqizPTCNLbKCd3oxWqzgR32YgFu8Jro7K0EL9zHZ+puYWmYBuf+ucNROOy8DMYC3Lrc7fij/p59IJHyTPnsa+lj6buEHN89owLWgYJ98PRl2QGDsBbgyhYwIKyEu7wXo4JPZ967rbBAsLn657nV3t/xVWLrmJd+boRp+oJRNnX3IvXYco6gwTIXZeOg5xiTPD2kiv4Z/1f2dbyBu3Bdp7Y/QRvq3obi/NTtxhv7AkSjiZysvAA1HjLKbfP5fn6Z9nasjXjrDMM3WTbc9IgJTubuiR2s3TcyMgxwV3OaR5p6/Xtrd8lFPdzfsV5Q88n4vIG2+xKWeg85qltRvLsRo52+lO7gFjcfHLhhwnHw/xk+w8BeGDjQxh1Ftavvnn84LB9n0wAFC8d1/1jOGaDnmVlboKROHua+kZcJ73ewFfP+xZmnZE7X7iDyIAO/cnDT/LbN3/LBcUf4D/mX5CevMvfDs1vSEeTklPGDeqNeh2rKvMw6ARbj3XTE4xy7al38N7yC/jxwd/zl51DOvAXG17kwdce5JzS8znN81/odYJVlXmZbbnHY1D/2tBuXtGSETdsFqOeU6oLiHjm857qDSyyl3Hni59nb+tQ4d7Pd/2cPx/8MzeuuJFLqi+hrS/M9vpubCYDy8vdmWXAkyQScu46thECnTJTX7AAT2ElHytaw1pLLQ9veXiwEUhPuIfP/ftzlDnKuOv0u0b+iAmN7XXdxDWNZZMdz3AiAeg6iqXrAJ9a+FH80T4eevVRAL7z+new6C18fNnHU54iHIvT0BWkyDX5QufjOav0HJoCR/jtvt/jMXtGWfNNRCgaJx7XsOdCdjZDSfXJuB7YAtQCWwf+vQX4f8B3c/HmQohLhRD7hBAHhBDZtZObKrL0eB6Ow2zIPHh2lXGGVy7+OzpeY2XB6ViNA2NJJGQhi9EmXSwyxGYyUOy2cKwzkNL2CZ2eq5Z+FI/Bxndee4gHN32bqBbhtlW3j28t1vzGmAU2E7Gw2Dm03XicF/UVyz/GJYVr+e6ex9nV9Bo/2/lzXmp8kf+s/CSXLVidXvfGjoOyfXledVqZrjKPlUWlLrr8EbYe68JhKeY753+b/niI2565iXDEz893/Zw/HfgT76m5hnLjWRS7LZlnSHobpTbcmgeVp4Np7EB3rs/B3CInVveFXF9xDa917uYrL9xBQkvwxZe+yJ6OPTy47kGqXXPZUd9DfWeQ6gJb5hnwJIFOKSGJhaB0FVSfBb4FUDAPQ9VprFzzLq4veC9H+o7ywIt30+Jv4YsvfZFaby23r7595KkiMbbXdw8GGVm1y434ZYa+/lUZfHUe4tOFS3HpHXz5la/KNu/xCDefcnPK08QTGofb/bhtxowLkMZE0yDUw+m+0zjYu4tgLDihL+tY+MMxLEb9pOymRpGDmg2QwXN8wLs1bfRG8gtqqbWVsqX1VQzCyEVzht1wdx6S4ytcNGGh81hU5dsJRxM0dKd2fplTfT7v8K3hf/b/lj/s+zOvtr7IO6s/wvyCcXYE+lqG3IGc6e8a5NlNzCt00NIb4kDryB20Qlc595+6gb3+er798r0c6z3Gl16+m2rHEj625Mb0PqPBbnmDbbJD2aoJEyYWo57VVQMB9NEu2vrDbDj3QdZ6FnLP1od5o/4lDnYf5I4X7mCOax5XlH4ai8HAmipvZtv/xwfO41xzUNr3AAAgAElEQVQzm8nAykoPmqWYq+fei1Nv5aZnbqa1t4Hn657n4S0Pc0n1Jdyw4gbqOgPsqO/GYTayqjIvbRvLkeOKynF1H5Nz/tzz5c1Q/lwoXkbxyktYP/86CvRuPv3MejqDHdz7yr20B9p5cN2D2I1DN9SJhMauxh76QjGWlrqzy6ompS2H/yXX8O6jnKULc7FrDX8+9Fv+fOB/efro0/zXkv8i35qf8lTJrHNNLpIRANEgby+TNRqvt21hbfHajDtVJuOckznzPO5PpmnaI8AjQohbNE17NNdvLITQI4Pwi4F64DUhxJ81Tdud6/fKisFim+wXH6fFMOAjq6UfOOh0zC87A+sOE8FEhHPLzxt6rm0PRANQvnbSTiDzCh209smJfmmZe9zj7AW1fKzqMr558HcItnJx+XtYXTZ/7IM7DkJ/C/gWpq1hTKLXCZaXe9h8tJNt9d0sLxtyqRA6HV9c9zW2/fk9fPCpawFYkbeOT6z8cHqBz0CDA5wlo1w/UlHmsWLUCXY19bLxUAfV+fO5d/UG7nztHj7+t6vY3nOQ1fnncab3Q1QX2Jnrs2cYODcNBc7la2CCdq1zCuyYDDp0ug9yub+e39c9yau/u5S6QBPrT7mduY61bDzUQSSWYEGRc/IZ51CPdCIxWOS4xsi+OTwF/MeZ17Pr2SP88ejfea19B+F4mAfXPTjCUL83FGV7XTcasLLCM7mFMEmgUwYQILNbzlIQgoKeOq5puYBHmv7Mvp6dvHv+u0e0zx2Lox1+wtEEy1L87adNqFfq1CN+LjW5+R2gE3rWFK+Z8KXH0xeK4bDkaOGJhcCS/c/nHBhPfyiWWYYrr5rT3fPZG2hkcd4qXOaBIDHQKecKV2nG80SSAoeZPLuJQ+1+it2W8f+ujBY+ueIG/vrMddy98S7yTIXcsuaasY8N98mbf7Mro3kiSVW+nVA0MRDUwIIix+B8cN78K/hA/b94/Mhf+UfzFoSmZ/3ye1laOk5jqeFE/NCwWSYjhmmcJ8Jmki3Ft9d1s6OuhyKXhbvPfoTrnr6S65+/DaEzoRcmPlx1D/l2J8vLM6xDiMfkuEI9ULpyTBnccFwWI2urvWyv03Ft9Rf47wN3ccNTH6ch3Mki7yI+tfKLbD3WTZc/Qr7DxLIy9+RuIuNRKSEJ90tpi2uMDXODiWUrL+K2cA93v/llzvvN+Who3L76dpYWDLVnjw84a7T1Sfmbz5nFjXY8Bo1bIdAhdz69c2V9QLCL2+Jx/rVlG194aQMes4erF6coZEXeZNd3BShxW3PgBx+XGfreBuYDRQYvLbHOjP2dk+OCHMnOEnEpO8sm2TIFpJJtXDDwzwYhxLuP/8rBe58KHNA07ZCmaRHg18xE/+gcZ54TCTJrlgLovHP5xtwbuc77Lq6oGsjadBwcuJueA/Y0vUDHwGLUU+mVVnHdgRRdB4Xg/afchEfvwKozc8faG8Y+bniAOoFf63iYDHK70WrUs62umzdb+gYz4w6rjzuWfQ6LMFFmKubeM+6m1JNGcNjbKCeGgQYHmX4QC10WTpvjxW01cqC1HyNnc4XvXWzrOUCVuZJrF36etdX5zCt0ZBY497dKXaXVk1bgnKTMY2VttZd3L7mTtfaV1AWaONN7AVX6y9jb1IfZoGftHO/kA+eIX2ZsdAZ5c5Zi29rrcvDp8+5nnrmCen8D1y28mXKH3AmJxBIcbvez+Yi06VpVOXmnD2AooNcbpbTFUyk9hXV6yKvmqjNuY5mlhhpzGdctSb3dGYjEONoRoNBlzrqIcvB6JWJQvIw18y7EpbOzwFqF05SZfCee0AhEYoPBatbEsq/ZAHCYDAgBvaEMd89MdlaVnAPAOu+AhCvUKwt2jVYoXJLVuBYUOYjGEqMyvcdTVnYqHy69lAK9m9sX34THOkamLuKH+s0yGZFGZjfVmJINVDYf7RrcRdM0jU+s+iJlpkKaQ83cuOgzXDB3wcTb/9GQ/PuCgc9jZr9Pi1HP2movc3x22vvD7G8xcHXNl4gnIBIL8fF5d7OmooY1VXmTC5yD3QOBc3pZervZwGk1+aydexZXV3yS/f56jJqB91d+id0NIfrDMRYWO6Xz0WQC56QFacQPZavHDpwH0OkEl576bq4u/zACWOtZw3vnyZ4FmqbR1hdm0+GOwcB50nMqDBSqb5U3jsXLoXiZ3GUUAmxeKhZdxnuKLkRD48Nz3o/DlHo3Yl9Ln+xMW5hl1lnT5I1/b4Nct0tWcN+cj3K99z+5tOyCiV9/HH2hGGajLrskSZKWXXDo+ezPk2NSzc7rgGeBd4zxnAb8Icv3LgOGt66rB0ap0oUQ1wHXAVRWpi46mxJMDsifn1V3riTJTFJ/OJZZAKHTYS+4mLdpW3A0bYZWo8yIO4pkdjdLqvNtNPeE2N3Yy6lzvONOVk0BC5+suZsa7TDFXQfB6h6axBMJaZfXcWAoQM0Ci1HPmqo89rf2c6wjwLGOAAa9IJ7Q0OvW8IPlP2SBrgFnoh0oTH2y4ZndLBZEm8nAKZV59IdjdPSHqfB+msUHiznX4qXcF4NMi838HTKAsLhkJinNwDmJ02Jk7RwfCwp/SNOef+DQ9MTydOR50uu6NS7RkLTwA6g4Na2F2ufx8J2Lvs9zW3/Jwj4P/9rXgsFgIDpgzehzmgc7Sk6aeExeL50BKk4bc1wWp4/vXfADjmx/nuDRelg+dhGtpmnsbuwFAQuKMtSmH0/SuQUNyk8DswM98LnaOzB0NpDobkDnSb9BSn84hqaRm+A5HgUtkZObf51OYDcb6EvVXGkcqoou5eO+Q1zpqoKDzw1ISUwysNFn93M6LUYq820c6whQ6DSPuwvVF4qypvQWrnCuY77VIKUZw7Ok/g4ZQKBB+alp6ZzHQwjBgiLZcW5fSx9bjnZJbbWAeFzjhoVfw9D3Apd7CyARAV2K308kIAPnWAQq1krJxiTQ6QRzfQ4qvTba+8NUhU7l8fwncHXvpLDQhSHflllSYXjgXLIiI3kLyF3GBUVOqvOvpXqbmaJgCJ9dYCt047WbJi9ZSgaCwS4pG3RMrKU36nXcct6nuWL3Gnpam9m69whYPSQ0jURCeu6vqsrLvsC5dbfMOBcvlzUtx6PTs/68u/E9n8cpwRJCoSAWy9h/h3WdATr7IywsznJeBSmf6m+VOy0DzlPLlr4L62tPEarbi2vR2RmdrjcUzcrycwSxcE7mr1yTatZKlvH/RNO0zJvBT8xYn9JRVR+apj0GPAawZs2azHt7ZovFNaFpf7okMzd9oRhFGZwyFI3THTPjm7sOdG0yk2T3ybvpHGxlGAY6/W052sWuxl6Wl4/Wozb1BDnc5mdF9RkscZ4is6WH/yUXH50B/G1S4+wqg6KlWTWUGT6uRSUuqvJttPaGCcXiGHQ6vHYTXnsRtNqlLlGLy+zV8e+paXJSaH9zIHBenXGAOhYOs2FQy7Ww6JMyw9EyoD1Pdxcg0ClfZ7QNbMFOPoBw2224l18Cx16G0B4wnznpcw1qBBMxGThnsFCXFVRx1ZmfwH94E526BvyuBZj0Ogqc5txMpC07ZWA/QUDv8ZWQX1FLR91ejhwtobqqetQxb7b00x2IsrjUlbW1E12H5VZ/6SowD2WKTlt4Gcdef4b+hl243Ol/VpPBqdOcg2uWgwYpw3GYDXSl2qEah45AnNMXfgJnfkAGD0a71KBm6QCSZK7PQXt/mJ2Nvaytzhu1hR2KxtlW143eaKJy+bmI1u0yA2j1yvk93DcwLpucJ8y5cVEodlvwOc209oXoCw3dFBU6V2MIz5VZ7rpN8j3H+qwFOmUgmEjIwNmahrxjAox6HSVuK7iBolroNEPbXpnhK1qS3t9pPCrnr2Tg7CqZ9HhMBh2Xr75K/j76j4G+EPST302lZZcMBAsXZTQunV7HnMVnE7W8SJf/IF2OU9EZTLitRnxOc3Y1GgB9zdLP3FszduA8gNVs5wNnXM+hzU/z5huvsmjVOaMyuN2BCPtbZYOp8rwsE3sRv0x6OYtHWLY6XR7ieXPp7zxMYaAzbWlVLJ4gEI5TPFlnp1EnDI1bAzSdpIpwkoKw1G7mk6ceGJ4WKgcap+i9ZgTJzE3Ktthj0NYnddcFeW7pl1y6cpQTQ7Z4bCbmF0qj/x31PYOduxIDBVW7GqRV0KJilwyYq8+Sj4EOKYlIBoEly3MSOA/HZjJQXWCnttjFvELH0N1/4SKpF+uph6MvygxzPCozIn0tsro6KSEpX5t1hmtMhJDZDbND6nD7Wyd+TV+LDFANZjmuXAQQJpsM3mIhmQVNZNCMJ0kiLhfzaECeazI6WUch9qJ5VIg2ah0hanyO3ATO3XUDDinz05rEK+avwON203F4O3ubhqzDEgmNvc291HUGqMy3UTrJZguDxCJSQuUoGqX1zHdYCLvm0NvXJ2sA0qQvFMOgF1hz4rSRfYOU4bgsRsLRRGp7y+MIx+L0BKL4XBYpsyk9RRad5ihwBpnFXFnhAWDL0a4RErQuf4TXjnQSS8jOmharHSpOl2478bD824qF5A5j9dk5C5yHj63EbWVBkZOFxU5KPVaZVbV5pVQrHpEWkO0HZJY5kZCyluadcgdI6OQNYw4C5zHxzpEBXU+dzIxO5KYSC8uAP9STdeA8iBBQvEJe+8bX5c3MZGjfPxSg5lVn/nq9EWPFKgptOhZyhPkD9p5ZB87RoPx9Wtzy72wCHC4vFXOXoPU2sn3foRFGA619IV6v68Zi0GfVYGrohHvl35hv0ainXCVz6Y1oxDpGdxgej+RY0yrgT4ccFTznmlTRxB4hxBHAJ4TYMez7AtA0TctuXx5eA+YLIeYADcAHgSuzPOeMx2Ux0toXyqhosK0/jM2kn3LD8aSWa39rHy8fbMdlNeIPxwhHExS5LCwpdQ1p80x2OXFON74FclFp3TWw5ToMg3n87bFcoh/QBddvlgF0Ya2s0j/+9zuYCd8/JNXIYQCBzSuz/s075CJYvHTi1yRJJOTYk9Xy9tQV3ikpWChvqpp3QNXZWTcYItwnK9NtBenr6HU6qmpXYdr9MofqDtDcW4nLaqQ/FCMSS1CVb2NeLjxau4/KnY+CBaOe0usEjvxSuo7sp7z7WNrb2n2hHOudIXeZ54Fx9YVimB3pBfft/TKQzarIKg1sJgOrKj3sqO9h85EuPDYj8YRGXyiG1aRnTZVnaEHX6aTjQobe0rkftBeqzpJ/3x375VcSoZNeyQULp+bGfzi+hXJ+6josd3eKl409NyUz4fGYvMFOQxKRNnqDzMAfe0XOpZVnZDZ3dB2VGVR3eXZyRotbfp7b9spz5mXuZDUCTRvwC9fG7bI7FnkVi9D1N3G4Yy8bNQduu4l4QqN/oJh4ZYUne7lGqBf8rfLnHeNa+1w29ljL6W9vwFO0KK2dyN6gDJ5d1hwVC+bAKngqSOW28SEhRDHwD+CKXL+xpmkxIcTNA+fXAz/VNG1Xrt9npuG2GWnsDhKIpNd9JxpP0B2IUJF3YrYtKvNteB0mjnb48YfjuCxGykqsWXdRmlIcPrCfKwO2cK+crMwusOXnPAs+Lgaz1OE2bZcLYW+TnHQtHkCTi07XYblF5iyRi1MOJCSjcJfJ9+g8KIvq0llEkprdQLsMvieolp8QnU5m44++LF0LytdMfpckEZd+4Tq93NXI4DzCWUxpaTnu3nYO26oJxDTybCbK8qw5ac5CIi4XV0fhuNlKn8vCYXMpge4WbEX+CRefeEKjPxzNusXuIDl0CwJwWaT0rCcYTXtOaOsLYzXpc5eJSoHTYuTUOV6OdQbo8kcw6KWuttRjyY3t31RgtMpajIhfejgnYvL3ZS84sVrPwloZQLXuhSP/GihG98l5JOKXGd2+ZjneytNy4uAyCqN1IIDeJPXU5aeml2DoaZAJA3uhnMOyJa9aridte2VyJhvpZseBAf31isz06jo97vLFLBHbaLT104YXo17HwmInZR5r9h7TIJM5OsO4Tcw8NiNxVzk93Q14+prTutnsDUUxG3XZB/aQU8OGXJMyetM0rRkYTC8KIfKACk3Tdoz/qvTRNO1vwN9yca7ZQrJ7VE8wmlbw3NYXJpFg8p3hJoHDbGBJ6RRMjFOJEHKxycJ5JGv0BihfLbeBOw8OZBuGYXKk9EDNGQXz5d165yFZLOarHT/ojMdkJsnfJvWOnvS7VKbE7JCLccsuqUsfpqXLiNY9EOkfkLdMYgL11WIPvMRSazcUjN9eeVL01MnrnDf+z1bgMLPPVkR3sB5bGotPbzBKIgFua452JGIhGfzk6CbSoNfhMBvoDqQnPYvFE3T6w5SfoJt/kJreuT4H5DApekIw2SddDJgz8qqlDrxtn5S8tb859JzOIGVy3pqpzYRb3PJmomGLlIdM5C7SeXggyPVKSWMu5IxCyATHkRflPF515uSSHYFOGTy7ylI6foyLqxRj5yGqtEaqqufm1q4tGpQ3Q3lV49oeCiHwelx0dNqp6GlAl07wHMxxsSDMrsxzEiHE88jMswHYBrQJIV7QNO1TUzy2kxK7SY9eL+gJRtPSW7b2hbEY9em1bFXMDDwVcusw1C09RoWQmfAcFZ5OiBAyEEbIwDXil/8/3j0g2CUzw5GAzNbkKnBO4qmUQXn7m3IXINOfv7dpSL842Zsii0tmhruOyMAgV4t+IiEXbWteSg22yaDD5XLS2WeltL9lwuC5OyiD0px93qdAL+ixmWjsDqYlPev0R0gkwDeTd64UI7G4ZHFiJCBlXImorGmx5k3NbtlY2Atk0NywRRZCFy0bLRGJRwd2+RoG2pSfktudRoNZZovrX5Pvk4kMDmQ9RNM2ee0KU3c5HRchpKSicav8OVN0xM2YriPycQJtuM9p5k2Tj/6+RlzhPjCP704UisYJROK5u1mewZnndP7S3Jqm9QLvBn6madpq4KKpHdbJixACt9VIT3DizE10IGtT5Jp5fziKCRBCLjbJQPpEBc7DKVosJ+1Ah3RGadouJ8zOQ1D3miyoTMTkIpXrwHlwDMtkxqppu5Q5pEvEL901LJ60CmxSkj9PBgDdR7M7z3D6m+XEnoYG2+ew0KPPJ9zfOeR+MQ49wSg2sz6ztsipiAVzYrM5HLd1QEucRrfU1r4wBr3Ak411omJ6MNlkMaCnUgazJypwTmLzSimcziAlHMc2yuLcrqOy+O7Q87JYPX+e1F9PhUTPXiB3lnrqZKF3uiR1zvGozIZnc9PuLJLJl46DExdzpks8KndIXSUTWjLm283E7EX0BqMTFsQnd6Q89lzd/Oe2ZiOXpPPXZhBClADvB/53isfzliDPZqI/FJuwYn06JBuKk4y8KpizTm4Z9rfJDErbPoj6B9wF1mVXHDgRBpPM3kT8MguTzuQfj8mMkxBy4cl2UbS4pXaz67A8dy7oPCRlOPaJtQGFLjMRS4FslJFi8dE0je5ABE+uJBsw0CAltzffyUC42586ARCLJ2jrC1OUC7cCxVsTi0sWHRcukq4k7W9KbXNfk/zsVZ0pZWpT+fdVsEAGr81vTHjzO0j7fllDUrg4N9rw/HnSBam3IftzgWywpsVTSs6S6HUCr9tJR8yC5m9LeWxXQNYZOHNlbhALg9Cn3U3zRJLOqnQfsqjvgKZprwkhaoD9E7xGkYJ8h1wcO/pT+6W29IaUZEORPUar1O/NuxDmDnzVnCc1wFNdxQ8ye+NbKAPHtn2pj00kBlpcD9jlZdGoYgT5c2W2pedY9ufqb5MOIN45aS3aFqMel8dDR1gvdwHGoTcYIxbXclPMCPJaxnNfqW4x6rGZ9HT4wymPa++PEE9oufN7Vbw10emktGDOOph3EdScLx9LV56YHT2dTr4XmnQBiU+wa9xTL2te3OW529HLZfY5kZA7kHZf2tfP5zTjN3gIdLenTEB0BSJ4bKbc3SxPwc1/rpgweNY07beapi3XNO3Ggf8f0jTtPVM/tJMXp9mAyaCj0z9+8ByKxun0Ryh2z8w/HMUsRAiZCc6lPV66eOfI7d+uw7Kaf6wFINlBMFm8mKYpf1pY86TVXefhzOQjY9F5SE7ozvQLgIrdVvp0LvzdreMufm39YYQYurnOminUCxY4zXQFIoP+2WPR1BPEYtQryYYid+iNsnjwRO9kmOzyZj4akEWM0dDYx3XXSUmJrSDrtvOjGMw+Z9kOo7deZvHTyDonKXCYiVu9skHSOAmAQCRGIBzHa8vlztnM9HiGFMGzEOLOgcdHhRDfPv7rxA3x5EMIQb7DRHt/GG2chbSpJ4SmkX0TB4ViplC4eCiArntVdicDGUz2t8nCIH/bwHFToMHOnysXjZ66yZ8j2A3BTrnwZCAnKXSaSVi9dPT6pZ3iGLT1hfHYTKO6iU2aZPCcq+z9MHwOs6yZHCcBEIwkb/6VZENxkmDPlzZ60aBsytVTP3QjHvFLS82WnXKnrWwKNNjOIlms13Fg8tlnTZMJBIs7I7meUa8jL7+YzmCMeP/Y0o2WXrkTVZjLGq0ZnHlO2SRl4HHziRjIWw2fw0xTd4j2/sio5gGaptHQFSTPbhrVZlahmLUkXUAsbqm9PvYK6IzSTk+LywxD+dqp02DbvNLOqvMQuCsnt7h1HpJjdmcW3Bv1OvIKSujev5uSvnaMx+kgA5EY/nCMMs/4lewZM4WZZ4/NiEEvaO0Lj1mT0dAdAMi+dbBCMZOwF0Dl6TK73PyGtOLU6aWUQ+hkdjh/3tRlxvPnDzhvNE6u+Vd/q8xeF6zM+KVleXb2G9x0dbRQUDLaeaS1N4TbZsRizGFh6QzOPKdqkvKXgcfHT9xw3joUOMyYDDoau4OjgueW3jChaJwFRbltE6tQzAjc5eAohr7GISs/i0e2uJ7qpjb586D+VZl9zrRzWMQv22x7505KK15emMebh2y0tTZT6htpWdfUIwPdnGdtYEoWHyEEhU4LLb0hYvHEiAYksXiChu4QPqc5twupQjETMDtlAB3oGGpqY7JL//4p2OUZgaNwKPvsKs08SO88KMc4iV4DeXYTRoeXtvZD5MciiGHyv0AkRl8oxvxcxiyxiEysZNuhdoqYcKUSQqwRQvxRCLFVCLEj+XUiBncyo9MJSj0W2vtloJxE0zQOtfXjsBimvJ2tQjFt6Ae6WhUtlpX0rpIT0w3Sni/1z52HZOFMJnQclNmlSbbrdVqM2N35tLU1E4kNvXc8IXea8h2m3GdtprBSvSzPSjyh0dg9Uv95pCNANJagyjvNDT8Uiqki2ZSrsFb6P3vnTH3gnHzfpPa5rymz1/rbpW+3t2bSmfHykhLCsQStbSNt+451BtDpoCiXxcEz2KYO0nPb+D/Az4D3AO8Y9qXIkjKPNBI/1OYf/F5dp2zdXeOzK62gQjEV5M+TE3Mmtk8Rv9wq9VRmJYOoKCuHWIRDjUOWdcc6A0RiCarzcxxsTrFe0G01kmc3caTDT3SgcDAYiXOs00+x24JbFQoqFLnHMUntc8cBOR+4Jt9opcBXhNWkp76hcfAzH4jEaOwOUuyy5vjmf+Y2SIH0guc2TdP+rGnaYU3Tjia/pnxkbwGsJj1V+TYau4McafdT1xlgf2sfPqeZQufMvNtSKGY99gIpE+nYn77vczLrnEZTlJRv7c6nwGmmrbWFQ239NHQHOdzeT6HLTF6uLOqSRENTng2bX+QgGk+wo76H9v4wr9d1oROCeYVKcqZQTAnJ7HPEn372OdApO8p6a7La4RN6I2XFRWjBTrbXddPRH2ZHfQ86IajxTcHNP+S8yVOuSEe4d7cQ4sfAM8CgsaemaX+YslG9hagpcOAPxznQ2g9Ant3IktJp6EanULyVKKyVHcs6D0oP6lSE+2XWOa8q+yyIyUFJnpNAMDi44+SyGllUMgWf+Vgot3Z/Y+CyyLHvaeqlyy8bJKwo9yits0IxlTiKZJOmjgOyfiRVQKxp0l/fYM640HnMt3YXUBM4wvZglNePdaPXC5aXuXP/mY+GADFjM8/pBM/XALWAEUgK9TRABc85QKcTrKjw0BOIoqHhthqVXEOhmGqsebLgpuuILGA0pciatO6WLYKzzDoDIAQ6q4f5pjBFxV5icY082xR85jXthFWql3qseO0m+sMx3FZj7qz2FArF2Aghb/obtkjrz/y54x/b2wihbtkoKxct1i1uvFYdZ5Rb6NcsuCxGTIYp+MzHQrInwQyNh9IJnldomrZsykfyFkfpAxWKE4yvVlo3Ne2AitPGzt701Muq+sLFucuAWNzQeQiXST91RZLxCKCdsGIbi1Gvss0KxYnEUSgz0B0Hxu8WGA1B2x4pU3NNwtpuLMzyfSxxPxZXDlqPj0csPGMlG5Ce5nmjEGLxlI9EoVAoTiQGMxQtlVmZtj2jnw/1SB9Xq1cWCuYKiwvQxm2WkhNmeKW6QqHIAUVLpJtO41Zp7TacRByatklXoeJlucvgmp3SxSfUk5vzjUcsOGMlG5Be8Hw2sE0IsW/Apu4NZVWnUChOClwlsltg9zHZ9CBZQOhvh7rX5MJUujK3W4cDmZspDZ6jM7tSXaFQ5ACDWbYNj4Vl06nQwJwSC0P9ZlkkWLwUzDks4BVCBtBTHjyHT4z93yRJR7Zx6ZSPQqFQKKaLwlrppNF5EPqaZcAcDcqCnLLVuQ9ATTbZpTCkMs8KhSJLrB6oOFXqn4++BEbb0Oe/ZIWs7cg1FreUtGna1GiS4zHZfGYG3/xPGDxPhS2dEOIbSK/oCHAQuEbTtO5cv49CoVCkhW8BOHyyuCYehbxqWZmeiwKbsbC4pli2EWYmV6orFIocYs2D6nWyc2q4FwxF4KlIXQidDRYXdMelXV4us9pJYkH5OINv/jPvMZsbngY+r2laTAjxdeDzwGenaSwKhUIhFyBr3ol5L7MLuo9KPeJUFA0m9YIztFJdoVDkGIMptetGLrEMFAqGeqYoeB5wRZ7BwfO0eAppmvaUpmnJ7gQbgcm3vFEoFIrZhsUFWgIi/TtP7oYAABMkSURBVFNz/hNkU6dQKN6CmByyaHCqds9mgexsJhhyXgv8fbwnhRDXCSE2CyE2t7W1ncBhKRQKxRQxPHMzFUxxa26FQvEWZrBocIqC5+hbOHgWQvxTCLFzjK93DjvmLiAG/J/xzqNp2mOapq3RNG2Nz+ebquEqFArFicNok41XpipzcwJacysUircwZieE+6bm3LGQLNyeKh/8HDBlmmdN0y5K9bwQ4mrgcuBCTdO0qRqHQqFQzDiSmZupWHziUdDiKvOsUCimDrNTFihGQ2DMcYY4FprRDVJgmmQbQohLkQWCV2iaFpiOMSgUCsW0YnJAeAo0z7NAL6hQKGY5Zqd8nIoEQCw04+ev6cqJfwdwAk8LIbYJIX4wTeNQKBSK6cHshER0SN+XK2ZBpbpCoZjlDAbPUyA9mwU1G9NiVadp2rzpeF+FQqGYMQzP3ORy2zM68z1SFQrFLEdvlHNMrjPPibiUns3wmo2Zq8ZWKBSKk5mpytwo2YZCoTgRmF25D54H56+ZnXlWwbNCoVBMB3qjXCByvfhEBxqkzOBKdYVCcRJgdsoug4lE7s45S2RnanZVKBSK6cLsyn2jlFlQqa5QKE4CzE5Ay+0cNktkZyp4VigUiuliKjI30eCM1wsqFIqTgGRr7lzuniWD5xk+h6ngWaFQKKYLs1O26Y76c3M+TZOZ5xm+8CgUipMAkwOELrfBcywIehPo9Lk75xSggmeFQqGYLkw5ztzEwjIYn+FbngqF4iRAiAG/+lxmnmfHzb8KnhUKhWK6MDkAkbtmKbHZseWpUChOEsxOiOQyeA7MivlLBc8KhUIxXeh0YLLnLnMTVTZ1CoXiBGJ2yh2vWCT7cyVlZ7Og4FkFzwqFQjGdmJ2583oezDzbcnM+hUKhSEUu23QnZWcq86xQKBSKlJidMtsSj2Z/rmgQdEbQT0vzWIVC8VYjl82eZpHsTAXPCoVCMZ3kMnMTDeW21bdCoVCkwmCWDZ9y4fU8i2RnKnhWKBSK6SQZPOdi8YkFZ4VeUKFQnETkqk13NCAfZ4HsTAXPCoVCMZ0YraAzqMyzQqGYnZid0jFI07I7Tyw0a2RnKnhWKBSK6cbszD54jscgEZ0VW54KheIkwuQALT6UOZ4s0eCsuflXwbNCoVBMN7kInmfRlqdCoTiJyFXdRjQ4K4oFQQXPCoVCMf2YnZCIycVjskRnT6W6QqE4iRgMnrOs24gGZ83NvwqeFQqFYrrJReYmmXk22bMfj0KhUKSLTi+D3mzs6mIRKf2YJbKzaQ2ehRCfEUJoQoiC6RyHQqFQTCumHAXPOqO0jVIoFIoTSbbSs1kmO5u24FkIUQFcDBybrjEoFArFjEBvkHKLbDI3kYCSbCgUiunB7JIBcCI+udcP7pyp4HkiHgbuBLL0NlEoFIqTgFxkbmbJwqNQKE4yspWeRVTmeUKEEFcADZqmbU/j2OuEEJuFEJvb2tpOwOgUCoViGjC75AIymcyNps2qYhuFQnGSYXbIx8kGz1G/7Fao0+duTFPIlDlRCyH+CRSP8dRdwAbgbemcR9O0x4DHANasWaOy1AqF4uTE7AQ0ufhYPZm9NhqUr1XBs0KhmA6MNhD6yXdKjQbBOHuKnacseNY07aKxvi+EWAbMAbYLIQDKga1CiFM1TWueqvEoFArFjGb4tmfGwfPs2vJUKBQnGULI7POkZRt+cBTmdkxTyAnvgahp2hvA4BUSQhwB1mia1n6ix6JQKBQzhmwyN7Os2EahUJyEmJ3Q35L56+IxiEdm1c2/8nlWKBSKmYAQky8ajPhl4D1LPFIVCsVJiNkJ8ShEQ5m9LuqXj7PIo37ag2dN06pV1lmhUCgY2PachF1duF8uPFIKp1AoFCces0s+ZpoAmGVOGzADgmeFQqFQDGB2DWRuMmzTHemfVVkbhUJxEpIMnkM9mb0u0g+IWTWHqeBZoVAoZgqW5OKTQfY5EYdYCEyOqRmTQqFQpIPeIAPg8CSCZ6N11tjUgQqeFQqFYuZgdgMis8xNssDQrIJnhUIxzZhdmWeew/2zbv5SwbNCoVDMFHQ6WXSTUfA8+4ptFArFSYrFDbGw/EqHRELOYbNs50wFzwqFQjGTsLgnETyLWdVgQKFQnKRY3PIx3TksGgA0FTwrFAqFIgssbkhEhzLKExHuG9ALqulcoVBMM+YM6zZmqexMzbYKhUIxk8g0cxPuHSo0VCgUiukkWTSY9vw1EDyrzLNCoVAoJo3JAUKX3uKTtLUzq+BZoVDMECzu9B03wj3S33kWOW2ACp4VCoViZqHTpV+xntwaTWarFQqFYrpJFg0mm5+kItQ7K+cvFTwrFArFTMPqkcFzIpH6uGQ3QrNz6sekUCgU6WD1ysdgV+rjYmHpUa+CZ4VCoVBkjdULWgJC3amPC/eCwSy/FAqFYiZgdoLOCMHO1MfN4p0zFTwrFArFTMOWZuYm1DPQWEWhUChmCEKANS+9+QtmZc2GCp4VCoVipqE3yuxNoGP8Y2JhaWdn9Zy4cSkUCkU6WD1yfkrVLCXUI5059IYTN64coYJnhUKhmInY8iHYPb7uOZnVSWapFQqFYqZgy5eP4yUANE3OYda8EzemHKKCZ4VCoZiJ2PJBi4+vGwx0gNAr2YZCoZh5WNygN0F/69jPh7plMyhbwYkdV45QwbNCoVDMRGz5Mjjuax77+UCnzNqozoIKhWKmIQTYfeBvH3v3zD+QkU5mqGcZatZVKBSKmYhODw4f9LfILc7hREOyra2SbCgUipmKo1Bml8dyDQq0y+y0wXTix5UDpi14FkLcIoTYJ4TYJYR4cLrGoVAoFDMWRxHEI6Or1vubh55XKBSKmYitQO6e9TaO/H4sIus5ZqlkA2BaShyFEOcD7wSWa5oWFkIUTsc4FAqFYkZjL5SLT0/9yCxzb6N04zA7pm9sCoVCkQq9AZxFUnrmqx1y1ehrBDRwFk/r8LJhujLPnwS+pmlaGEDTtHEU5QqFQvEWRm8Adxn0NUE0KL8X7JIWT+7y6R2bQqFQTISnUko3eurk/zUNuo5Kb2fL7PN3TjJdwfMC4BwhxCYhxAtCiLXjHSiEuE4IsVkIsbmtre0EDlGhUChmAN4a+diyGxJxaN0jq9hdKnhWKBQzHGuelGd0HIBIALqOQDQA+XOne2RZMWWyDSHEP4GxcvJ3DbxvHnA6sBb4jRCiRtOOr4oBTdMeAx4DWLNmzajnFQqF4qTGaIWCBdC2Fw78U7btLlk5KxsLKBSKtyBFi+HoK3Dk33L+chTOaskGTGHwrGnaReM9J4T4JPCHgWD5VSFEAigAVGpZoVAojsc7Bwxm6e3sKJKLj0KhUMwGTHaoPB26j4LBAnlzpntEWTNdqYs/ARcAzwshFgAmoH2axqJQKBQzH1ep/FIoFIrZhtkBRUumexQ5Y7qC558CPxVC7AQiwNVjSTYUCoVCoVAoFIqZxLQEz5qmRYCrpuO9FQqFQqFQKBSKyaI6DCoUCoVCoVAoFGmigmeFQqFQKBQKhSJNVPCsUCgUCoVCoVCkiQqeFQqFQqFQKBSKNFHBs0KhUCgUCoVCkSYqeFYoFAqFQqFQKNJEzCZ7ZSFEG3B0Gt66ANXEJRPU9cocdc0yQ12vzFDXKzPU9coMdb0yQ12vzJjO61WlaZrv+G/OquB5uhBCbNY0bc10j2O2oK5X5qhrlhnqemWGul6Zoa5XZqjrlRnqemXGTLxeSrahUCgUCoVCoVCkiQqeFQqFQqFQKBSKNFHBc3o8Nt0DmGWo65U56pplhrpemaGuV2ao65UZ6nplhrpemTHjrpfSPCsUCoVCoVAoFGmiMs8KhUKhUCgUCkWaqOBZoVAoFAqFQqFIExU8DyCEeJ8QYpcQIiGEGNcSRQhxqRBinxDigBDic8O+P0cIsUkIsV8I8T9CCNOJGfn0IITwCiGeHvh5nxZC5I1xzPlCiG3/v727j7GjKuM4/v3Zpm2IsbQ1apFi26QW8a2tTVPUAGIDSmJbY6NLQmwFYwoG/jAkQmqM0SjCP02MRkSDoCSgXUXXF2yAtjEi9QWxL1DabosJpZWKvMUQV6yPf8xZPd6du3dWdu/L3N8nubkz55yZnvPkmenZuWfvZq+/S1qf6m6T9HhWt6z9o2ifKvFK7U5lMRnKyp1fY9ssk/Rgum73SvpIVtcX+dXsfpTVz0z5MpzyZ2FWd30qPyjp4nb2u1MqxOtTkh5N+XS/pDdkdaXXZp1ViNcmSX/J4vLxrG5jun4PS9rY3p53ToWYbc3idUjSc1ldX+WYpFslnZS0v0m9JH0lxXKvpBVZXWfzKyL8KtZ9vwlYCuwCVjZpMw04AiwGZgB7gHNS3feBgbR9M3Blp8c0xfG6CbgubV8H3Nii/VzgGeC0tH8bsKHT4+i2eAF/a1Lu/Brb5o3AkrR9BnACOD3t1z6/xrsfZW2uAm5O2wPA99L2Oan9TGBROs+0To+pC+L1nuwedeVovNJ+6bVZ11fFeG0Cvlpy7FzgaHqfk7bndHpM3RCzhvZXA7dm+/2WY+cBK4D9TeovAe4BBKwGfpPKO55ffvKcRMSBiDjYotkqYDgijkbEP4C7gHWSBFwIDKZ2twPrp663XWEdxTih2ng3APdExItT2qvuNdF4/Yfzq3y8EXEoIg6n7ePASWDMX4KqsdL7UUObPI6DwHtTPq0D7oqIkYh4HBhO56uzlvGKiJ3ZPWo3cGab+9hNquRXMxcD90bEMxHxLHAv8L4p6mc3mWjMLgXubEvPulBE/JLioVoz64DvRGE3cLqk+XRBfnnyPDGvB57I9o+lsnnAcxHxz4byOnttRJwASO+vadF+gLE3iS+mj2K2Spo5FZ3sIlXjNUvS7yXtHl3igvOrZX5JWkXxpOdIVlz3/Gp2Pyptk/LneYp8qnJs3Ux0zFdQPPUaVXZt1lnVeH0oXWeDkhZM8Ni6qTzutCRoEbAjK+63HGulWTw7nl/T2/mPdZqk+4DXlVRtiYgfVzlFSVmMU97TxovXBM8zH3grsD0rvh74M8WE5xbg08Dn/7+edodJitdZEXFc0mJgh6R9wAsl7Zxf/z3PfOC7wMaI+Fcqrl1+lahy3+mre1YLlccs6TJgJXB+Vjzm2oyII2XH10SVeP0EuDMiRiRtpviU48KKx9bRRMY9AAxGxKmsrN9yrJWuvX/11eQ5Ita8zFMcAxZk+2cCx4GnKT5OmJ6e7oyW97Tx4iXpKUnzI+JEmrycHOdUHwbujoiXsnOfSJsjkr4NXDspne6gyYhXWn5ARByVtAtYDvwA51dpvCS9CvgZ8Jn0sd7ouWuXXyWa3Y/K2hyTNB2YTfExaZVj66bSmCWtofgB7vyIGBktb3Jt1nli0zJeEfHXbPebwI3ZsRc0HLtr0nvYfSZyXQ0An8wL+jDHWmkWz47nl5dtTMzvgCUqvvlgBkXyD0Wxgn0nxbpegI1AlSfZvWyIYpzQerxj1nWlCdHoet71QOlv29ZIy3hJmjO6vEDSq4F3AY86v5rGawZwN8WauG0Ndf2QX6X3o4Y2eRw3ADtSPg0BAyq+jWMRsAT4bZv63Skt4yVpOfANYG1EnMzKS6/NtvW8M6rEa362uxY4kLa3AxeluM0BLuJ/P3msqyrXJJKWUvyi24NZWT/mWCtDwEfTt26sBp5PD0Y6n1/t/O3Ebn4BH6T4aWYEeArYnsrPAH6etbsEOETx0+CWrHwxxX8+w8A2YGanxzTF8ZoH3A8cTu9zU/lK4FtZu4XAk8ArGo7fAeyjmNTcAbyy02PqdLyAd6aY7EnvVzi/xo3XZcBLwB+z17J+yq+y+xHF8pS1aXtWypfhlD+Ls2O3pOMOAu/v9Fi6JF73pfv/aD4NpfKm12adXxXidQPwSIrLTuDs7NjLU94NAx/r9Fi6JWZp/3PAlxuO67sco3iodiLdx49R/J7BZmBzqhfwtRTLfWTfhNbp/PKf5zYzMzMzq8jLNszMzMzMKvLk2czMzMysIk+ezczMzMwq8uTZzMzMzKwiT57NzMzMzCry5NnMrMdJ+lP6btiX1cbMzFrz5NnMzMzMrCJPns3MeoikH0l6SNIjkj7RULdQ0mOSbpe0V9KgpNOyJldL+oOkfZLOTseskvRrSQ+n96VtHZCZWY/x5NnMrLdcHhHvoPhri9dImtdQvxS4JSLeBrwAXJXVPR0RK4CvA9emsseA8yJiOfBZ4EtT2nszsx7nybOZWW+5RtIeYDewAFjSUP9ERDyQtu8A3p3V/TC9PwQsTNuzgW2S9gNbgTdPRafNzOrCk2czsx4h6QJgDXBuRLwdeBiY1dAsxtkfSe+ngOlp+wvAzoh4C/CBkvOZmVnGk2czs94xG3g2Il5Ma5ZXl7Q5S9K5aftS4FcVzvlk2t40Kb00M6sxT57NzHrHL4DpkvZSPDHeXdLmALAxtZlLsb55PDcBN0h6AJg2mZ01M6sjRTR+wmdmZr1I0kLgp2kJhpmZTQE/eTYzMzMzq8hPns3MzMzMKvKTZzMzMzOzijx5NjMzMzOryJNnMzMzM7OKPHk2MzMzM6vIk2czMzMzs4r+DYDJAcvtjlhEAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Frequency misfit term alpha = (f - f')/fs\n",
    "\n",
    "N = 10\n",
    "\n",
    "def misfit_term(alpha, N):\n",
    "    n = np.arange(0, N-1)\n",
    "    \n",
    "    alpha_n = np.outer(n, alpha)\n",
    "    return np.sum(np.exp( 2j*np.pi* alpha_n), axis=0)\n",
    "\n",
    "def almost_sinc(\n",
    "    x, phase=np.pi/4):\n",
    "    return np.sin(2*np.pi*alpha*n + phase) / (np.sqrt(2)*np.pi*n)\n",
    "\n",
    "def almost_sinc(alpha, N, phase=np.pi/4):\n",
    "    n = np.arange(0, N-1)\n",
    "    \n",
    "    alpha_n = np.outer(n, alpha)\n",
    "    return np.sum(almost_sinc_per_n, axis=0)\n",
    "\n",
    "####\n",
    "alphas = np.linspace(-1, 1, 2*500)\n",
    "misfit = misfit_term(alphas, N)\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(12,4))\n",
    "ax.set_xlabel(\"alpha\")\n",
    "ax.set_ylabel(\"misfit sum term\")\n",
    "ax.axhline(N-1)\n",
    "\n",
    "ax.plot(alphas, np.real(misfit), alpha=0.3, label=\"Real\")\n",
    "ax.plot(alphas, np.imag(misfit), alpha=0.3, label=\"Imag\")\n",
    "ax.plot(alphas, np.abs(misfit), label=\"Abs\")\n",
    "#ax.plot(alphas, np.abs(misfit)**2, label=\"Power\")\n",
    "if True:\n",
    "    ax.plot(alphas, (N-1)*almost_sinc(np.sinc(2*np.pi*1.5*alphas)) + (N-1)*np.abs(np.sinc(2*np.pi*1.5*(alphas -1 ))))\n",
    "if True:\n",
    "    n= 10\n",
    "    ax.plot(alphas, -1*(N-1)*np.sin(2*np.pi*alphas*n - np.pi/4))\n",
    "\n",
    "fig.legend(loc='center right')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Header"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "s\n"
   ]
  },
  {
   "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}