{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Beacon Sync\n",
    "\n",
    "Synchronise two delta peaks, by using an intermediate beacon that was sent out together with it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.signal as signal\n",
    "import scipy.fft as ft"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Monkey patch correlation_lags\n",
    "if not hasattr(signal, 'correlation_lags'):\n",
    "    def correlation_lags(in1_len, in2_len, mode='full'):\n",
    "        r\"\"\"\n",
    "        Calculates the lag / displacement indices array for 1D cross-correlation.\n",
    "        Parameters\n",
    "        ----------\n",
    "        in1_size : int\n",
    "            First input size.\n",
    "        in2_size : int\n",
    "            Second input size.\n",
    "        mode : str {'full', 'valid', 'same'}, optional\n",
    "            A string indicating the size of the output.\n",
    "            See the documentation `correlate` for more information.\n",
    "        See Also\n",
    "        --------\n",
    "        correlate : Compute the N-dimensional cross-correlation.\n",
    "        Returns\n",
    "        -------\n",
    "        lags : array\n",
    "            Returns an array containing cross-correlation lag/displacement indices.\n",
    "            Indices can be indexed with the np.argmax of the correlation to return\n",
    "            the lag/displacement.\n",
    "        Notes\n",
    "        -----\n",
    "        Cross-correlation for continuous functions :math:`f` and :math:`g` is\n",
    "        defined as:\n",
    "        .. math::\n",
    "            \\left ( f\\star g \\right )\\left ( \\tau \\right )\n",
    "            \\triangleq \\int_{t_0}^{t_0 +T}\n",
    "            \\overline{f\\left ( t \\right )}g\\left ( t+\\tau \\right )dt\n",
    "        Where :math:`\\tau` is defined as the displacement, also known as the lag.\n",
    "        Cross correlation for discrete functions :math:`f` and :math:`g` is\n",
    "        defined as:\n",
    "        .. math::\n",
    "            \\left ( f\\star g \\right )\\left [ n \\right ]\n",
    "            \\triangleq \\sum_{-\\infty}^{\\infty}\n",
    "            \\overline{f\\left [ m \\right ]}g\\left [ m+n \\right ]\n",
    "        Where :math:`n` is the lag.\n",
    "        Examples\n",
    "        --------\n",
    "        Cross-correlation of a signal with its time-delayed self.\n",
    "        >>> from scipy import signal\n",
    "        >>> from numpy.random import default_rng\n",
    "        >>> rng = default_rng()\n",
    "        >>> x = rng.standard_normal(1000)\n",
    "        >>> y = np.concatenate([rng.standard_normal(100), x])\n",
    "        >>> correlation = signal.correlate(x, y, mode=\"full\")\n",
    "        >>> lags = signal.correlation_lags(x.size, y.size, mode=\"full\")\n",
    "        >>> lag = lags[np.argmax(correlation)]\n",
    "        \"\"\"\n",
    "\n",
    "        # calculate lag ranges in different modes of operation\n",
    "        if mode == \"full\":\n",
    "            # the output is the full discrete linear convolution\n",
    "            # of the inputs. (Default)\n",
    "            lags = np.arange(-in2_len + 1, in1_len)\n",
    "        elif mode == \"same\":\n",
    "            # the output is the same size as `in1`, centered\n",
    "            # with respect to the 'full' output.\n",
    "            # calculate the full output\n",
    "            lags = np.arange(-in2_len + 1, in1_len)\n",
    "            # determine the midpoint in the full output\n",
    "            mid = lags.size // 2\n",
    "            # determine lag_bound to be used with respect\n",
    "            # to the midpoint\n",
    "            lag_bound = in1_len // 2\n",
    "            # calculate lag ranges for even and odd scenarios\n",
    "            if in1_len % 2 == 0:\n",
    "                lags = lags[(mid-lag_bound):(mid+lag_bound)]\n",
    "            else:\n",
    "                lags = lags[(mid-lag_bound):(mid+lag_bound)+1]\n",
    "        elif mode == \"valid\":\n",
    "            # the output consists only of those elements that do not\n",
    "            # rely on the zero-padding. In 'valid' mode, either `in1` or `in2`\n",
    "            # must be at least as large as the other in every dimension.\n",
    "\n",
    "            # the lag_bound will be either negative or positive\n",
    "            # this let's us infer how to present the lag range\n",
    "            lag_bound = in1_len - in2_len\n",
    "            if lag_bound >= 0:\n",
    "                lags = np.arange(lag_bound + 1)\n",
    "            else:\n",
    "                lags = np.arange(lag_bound, 1)\n",
    "        return lags\n",
    "\n",
    "    signal.correlation_lags = correlation_lags"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "### signal generation\n",
    "def fft_bandpass(signal, band, samplerate):\n",
    "    \"\"\"\n",
    "    Simple bandpassing function employing a FFT.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    signal : arraylike\n",
    "    band : tuple(low, high)\n",
    "        Frequencies for bandpassing\n",
    "    samplerate : float\n",
    "    \"\"\"\n",
    "    signal = np.asarray(signal)\n",
    "\n",
    "    fft = ft.rfft(signal)\n",
    "    freqs = ft.rfftfreq(signal.size, 1/samplerate)\n",
    "    fft[(freqs < band[0]) | (freqs > band[1])] = 0\n",
    "    \n",
    "    return ft.irfft(fft, signal.size), (fft, freqs)\n",
    "\n",
    "def deltapeak(timelength=1e3, samplerate=1, offset=None, peaklength=1):\n",
    "    N_samples = int(timelength * samplerate)\n",
    "    if offset is None:\n",
    "        offset = (np.random.random(1)*N_samples).astype(int) % N_samples\n",
    "    elif isinstance(offset, (tuple, list)):\n",
    "        offset_min = offset[0]\n",
    "        offset_max = offset[-1]\n",
    "        \n",
    "        offset = (np.random.random(1)*(offset_max - offset_min)+offset_min).astype(int) % N_samples\n",
    "    \n",
    "    position = (offset + np.arange(0, peaklength)).astype(int) % N_samples\n",
    "    \n",
    "    signal = np.zeros(N_samples)\n",
    "    signal[position] = 1\n",
    "    \n",
    "    return signal, position\n",
    "\n",
    "def sin_delay(f, t, t_delay=0, phase=0):\n",
    "    return np.cos(2*np.pi*f * (t - t_delay) + phase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Beacon initial [ns]: 897.9591836734695\n",
      "Beacon initial [idx]: 4489.795918367347\n",
      "Beacon difference [ns]: 35.71428571428571\n",
      "Impulse offsets [ns]: [ 65.6 188.4]\n",
      "Time difference peaks [ns]: 122.8\n"
     ]
    }
   ],
   "source": [
    "us = 1e3 # ns\n",
    "ns = 1/us # us\n",
    "\n",
    "\n",
    "band = (30, 80) # MHz\n",
    "samplerate = 5000 # MHz\n",
    "timelength = 0.3 # us\n",
    "\n",
    "time = np.arange(0, timelength, 1/samplerate)\n",
    "\n",
    "# generate beacons\n",
    "if True: # in-band\n",
    "    f_beacon = 70 # MHz\n",
    "else: # under band\n",
    "    f_beacon = 20 # MHz\n",
    "\n",
    "beacon_time_offset = 2.5/f_beacon # us\n",
    "beacon_amplitude = 0.1\n",
    "beacon_init_phase = 2*np.pi*  4.4/ns/f_beacon\n",
    "\n",
    "beacons = np.array([\n",
    "    beacon_amplitude * sin_delay(f_beacon, time, phase=beacon_init_phase, t_delay=0),\n",
    "    beacon_amplitude * sin_delay(f_beacon, time, phase=beacon_init_phase, t_delay=beacon_time_offset)\n",
    "])\n",
    "\n",
    "# generate impulses\n",
    "impulses = []\n",
    "impulses_offsets = []\n",
    "impulses_def_offsets = [(0.05*samplerate,0.2*samplerate)] # random offsets in interval\n",
    "for i in range(2):\n",
    "    offset = None\n",
    "    if impulses_def_offsets:\n",
    "        if len(impulses_def_offsets) == 1:\n",
    "            offset = impulses_def_offsets[0]\n",
    "        else:\n",
    "            offset = impulses_def_offsets[i]\n",
    "    orig_imp, imp_offset = deltapeak(timelength, samplerate, offset=offset, peaklength=1)\n",
    "\n",
    "    ## Bandpass it\n",
    "    imp, _ = fft_bandpass(orig_imp, band, samplerate)\n",
    "    imp /= np.max(imp)\n",
    "    \n",
    "    impulses.append(imp)\n",
    "    impulses_offsets.append(imp_offset/samplerate)\n",
    "\n",
    "impulses = np.asarray(impulses)\n",
    "impulses_offsets = np.asarray(impulses_offsets)\n",
    "print(\"Beacon initial [ns]:\", beacon_init_phase/(2*np.pi*f_beacon) /ns)\n",
    "print(\"Beacon initial [idx]:\", beacon_init_phase/(2*np.pi*f_beacon)*samplerate)\n",
    "print(\"Beacon difference [ns]:\", beacon_time_offset/ns)\n",
    "print(\"Impulse offsets [ns]:\", impulses_offsets[:,0]/ns)\n",
    "print(\"Time difference peaks [ns]: {}\".format( (impulses_offsets[1,0]-impulses_offsets[0,0])/ns ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "full_signals = impulses + beacons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAscAAAEHCAYAAABC9usHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOydd3hUVfrHP2f6TGbSG6GF3pEqTUXEtnZUsKzYsWDvfe2r6+qK4vqzL3Z3sSuKIoIiSAnSW0BqCGmTyWR6Pb8/ZhISSANSQM/neSKPM/fO97z3lPvec97zXiGlRKFQKBQKhUKhUICmrQugUCgUCoVCoVAcLijnWKFQKBQKhUKhiKOcY4VCoVAoFAqFIo5yjhUKhUKhUCgUijjKOVYoFAqFQqFQKOIo51ihUCgUCoVCoYija6kfFkK8BZwBlEgp+zflnPT0dJmbm9tSRWozohs2AKDp06fVNP1+PwAmk0lpHiibNsX+7dWr9TQbQWkqTaWpNJWm0lSazcfy5cvLpJQZdX3XYs4xMAN4CXinqSfk5uaSl5fXYgVqK7wjRgBgWbKk1TTz8/MB6Nmzp9I8UI4/Pvbv/Pmtp9kISlNpKk2lqTSVptJsPoQQO+r7rsXCKqSUPwPlLfX7CkVLIKXE4Q1RVOnH6Q21dXEUCoVCoVC0Mm0ecyyEuEYIkSeEyCstLW3r4ij+5Lz281Y2FVWyvczDhJcXUulXDrJCoVAoFH8mREu+PloIkQt83dSY42HDhsk/YliF4sggv9jFX15YwLefPUSmzcTg8fdz9TFdeOD0vm1dNIVCoVAoFM2IEGK5lHJYXd+1+czxnwH3Sy/hfumlVtW02+3Y7XaleQC8Mv93jDoNOSk6zEbJxKEdmLFoO2XuQItpNhWlqTSVptJUmkpTabYOyjluBTTvvovm3XdbVfPP0uibS7PcE+TLVbvp1fdHNpSvZXXpKlJyfiIUkXyyvKBFNA8Epak0labSVJpKU2m2Di3mHAshPgR+BXoJIQqEEFe1lJZCcajMXlsEtmVsCXxLhiWTVHMaH2x+nd5dCvjkt4LGf0ChUCgUCsUfghZL5SalvKilfluhaG6+WrONhOw5HJU5mM6JO5EySrekVBz6z9mxbSo77B46pyW0dTEVCoVCoVC0MPXOHAshdEKIa4UQs4UQq4UQq4QQ3wohrhNC6FuzkApFS+INhllRPpeoppKbB9+MADRCw/WDrqc8VIDOto4564vbupgKhUKhUChagXqzVcTDIiqAt4GqdeUOwGVAqpTyguYujMpWoWgLfsov5fofJ9M5zcS3Ez+vfglI5Me5nPrpqTgcafTV3s47Vx7dtgVVKBQKhULRLBxstoohUsrrpZSLpZQF8b/FUsrrgcEtU9Q/JipbxeGt+c3GlWjNu5nU+9xan2s1Ws7seiYBw3qW79pOOBJtNs0DRWkqTaWpNJWm0lSarUNDzrFDCDFRCFF9jBBCI4S4AHC0fNH+OKhsFYe35i+FPwFwerdT9/vuL13+AkiCxnWsLaxsNs0DRWkqTaWpNJWm0lSarUNDzvGFwPlAsRAiXwiRDxQB58a/UyiOeCr9IRxyJen6bmRaMvf7vntyd7Is7dBaN7Jk6+HfoRUKhUKhUBwa9TrHUsrtUsoLpJQZwChgtJQyM/7ZttYrokLRcizZsRONeSfDM8fU+b0QgnEdx6K3buHXbWpTnkKhUCgUf3SalOdYSmmXUpa1dGEUitbmu99/QgjJWT1PrPeY4zocByLIipI8WvJ16wqFQqFQKNqeerNVtAUqW4WitTn5vRvYE8pj1eW/oqkKr49nq2D+fAD8YT8jPxiNr2wUP1z+TzqmWtqkrAqFQqFQKJqHg81WoWgmVLaKw1ezJLieFNFrr2NcByadiR5J/dBatrK6wHlE2qk0labSVJpKU2n+GTUPhiY5x0KIFCHE0UKI46r+WrpgfyRUtorDU3N7RSERbRk9kgY1euwxHY5GY9rN8p2FR5ydSlNpKk2lqTSV5p9V82Bo1DkWQlwN/Ax8Bzwa//eRli2WQtHyfJ2/AIAxHRp/uceInOEIIVla9FtLF0uhUCgUCkUb0pSZ41uA4cAOKeU4Yi8AKW3RUikUrcDiwmXIiIkTuzc+c3xUxlEItOzwrCEaPXzi9BUKhUKhUDQvTXGO/VJKP4AQwiil3Aj0atliKRQtz1bXejTBznRMTmj0WLPOTI65B2HDNnZXeFuhdAqFQqFQKNqCRrNVCCE+A64AbgVOIPZ2PL2U8rTmLozKVqFoLbwhLyPeH0mWPIMfrvh77S/3yVZRxb3zH+frbZ/x6FGfcd6Qzq1SToVCoVAoFM3PIWWrkFJOkFJWSCkfAR4C3gTOad4i/rFR2SoOP82VJWtASPqk9G/yOcd2HIrQhPhh7cIjxk6lqTSVptJUmkrzz6x5MDQ1W4VWCJEDbANWAtktWqo/GCpbxeGnOX/7MgCO6TikyecMzorFJq/eufyIsVNpKk2lqTSVptL8M2seDLrGDhBC3AQ8DBQD0fjHEhjYguVSKFqU34pXEw2mcXTnjk0+p11CO4wiidLgVvWmPIVCoVAo/qA06hwTy1bRS0p5+Lv6CkUTkFKyw70BEehK5wN4250Qgk4JfVitX02R0692pSoUCoVC8QekKWEVuwBnSxdEoWgtir3F+GUF2caeaDTigM4dnHkUGoOD1Xv2tFDpFAqFQqFQtCVNyVbxJrHUbbOAQNXnUsp/NXdhVLYKRWvw7bbvuPvnOzkh8UlemHDW/gfUk60CYGHBYq6bO4XxKfcz7ayLWrScCoVCoVAoWoaGslU0JaxiZ/zPEP9THCBVmSqsN97YappVAe9paWnN+rvRqGTZ9nJ+21mBEDC0cwrDOqcghGgxzYY4GM2Fu5YjozpGdxhwwHqDswYSdkVY6VoCtJ5z3BzXtrjSz48bSyipDNApzcwJvbJIsuhbVPNAUZpKU2kqTaWpNNuaRp1jKeWjrVGQtsYZcPL5ls9Zb1+PWWdmZM5Ixncaj15Tv/PQVKozVRzhzvHa3U4e/HwtK3dV1Pp8eG4Kf58wAOk8MjraypI1RAPtOKrjgZfTordg9GZQ5F9/wOceCgdbn4FIgJXFa3hv2Xq+X+0l4MkBtADYjDpuObEHV47pUmd4yZ9l4KxPc6tzK78U/IIr5KJbUjfGdhyLWWduUc2WRGkqTaWpNP+MmgdDU7JV9ATuBHJrHi+lPKHlitW6fPX7Vzy19ClcQRc5CTl4w14+2fwJfVL78MQxT9AzpWdbF7HN+WpVIXfOXEWiWc/T5w7gtIHtkDL2+b/m5DPh5UXcMyqREV0P7wYfjoYp8GwG/zC6Z1gP6jfaW7qRH15GcaWPrMTmcZaaG2fAyWurX2Pmppn4Ij4ADJ0gTZ/EpX0vZUjy2bw8bwdPzNrA0m3lvHDhYMwGbRuX+vDAFXTxj6X/4Mvfv0SyN+ws05zJ/SPuZ3zn8W1YOsWRijfkxaA1oNM0ZcFWoVC0JU3ppTOBV4A3gEjLFqd1kVIy7bdpvLX2LYZmDeX+zmfRc/HrRIt2MCcxiad0O5j8zWSmjZvGqJxRLVOIvBmw+GXwO8GcDCOuh2GXt4zWvpqeUtBoYdBf4aT6Fwhmrd7DLR+tYFjnVP7vkiGkWY04/A7Wla8jKcPBPecJXp/r5qHPd/D3cwfSs65niV1L4YeHoWgdmJPgmDvqtbPCGyQYjpJhMyLEgW2Yq9POGtf2966jiBAg29QTnbZJab73Y0B6HzZ7fmHe1rVcOGh47S8PwM768ATCOH0hshNNTdswuI+dy/qfwb0l87H77JhDQwkV9+Hm44fTPSfAF79/wb9XTadHymyenfAsx63N4IlZ67n6nWW8edlwTPomOsi7lsLCabArD4g22oaahV1L8f3wEN9WbuE3swWy+jG4z0T+0uUvWPRNzzrSEIXuQq7/4Xp2Vu7k8v6Xc0nQQOryt1kecvCcrYhb59/KDYNu4NqB1x5a22yIZmhDB0xbjENtYWcraxa6C3l3yT/5vmA+JYTRS8kwUxaXHvsoY3LG/CHaUKm3lO93fM+Wzd8QLV5Hrt/LCSFJ5wEXt8qY8KfpKwun4SrIo0gjMfc5m5xTn0EjDu4e1mTNVri2s7fNZlTOKJKMSc3+24dCUzbkLZdSDm2NwrTFhryn5t5KuHAFdxXswBhw1PquRKvl+uwsdup1vDr0HoYOuOSgNLwjRgBgWbKEUCTEosJFrF/2b+Tu5fT3uhnj81PLJdFZIKXTIXW6/Px8AHpWeapVzsz2ReB37H+C1gAdhsGJj0LHo6s/nrexhGvezWNQx2RevXQgPxbM5uP8j1lnX7ffTwQLBRFPL56YdBMXDDg2NvDnzYB5fwdP8f6aBht0HQtjboGOR7Ou0MnjX69n8dZyAHpn23hywgCGdk5pup0Acx6Gpa9DyLPf8TOTUngs1cZfNddx7+Qb6v7RBjbkAXy56DvuXHATp/S9g+lnXtv4td3HzvpwB8I8OWs9nyzfTTASezi4bmw3rhidy5Ytmxu1Mwy8kpzEa8mJdAxHudOu5fvKkxh70d2c2n/ve3sWFCzgwYUPEoqGeGHcC+wqbMcdM1dxTPd03rhsGEZdrDU2Vxs6EOqsz3gbWhF2cG9GOoV6Henh2HN6mU5LZlTwWJdzGXP8I4ekmdkpk8nfTsbut/O0qT+j136FNuKtPi4MPJKeypc2KzdFE7nmxOdbxM6m9JXm0tz6ze2s3vQZpYQRQGYkQv9AgC6hMMJgg6T2rTsO6czQfnDzt6E5D0PefyAQS7zk1AjWGYxUaDXYhJ4B5mySR0xtNjullHz4xWSmOVYSEjDO66NPIIhDq2GuxUKhXsfJ0sxjY/9JQpexzWdnvA0VBMrYodejl5JewRBJ0Wiz31e8IS/TV0znow3vE0aSGokgJNjj48dJHi+3VXjo2K4F6rOhvtKSbaie+0pL3rPDu/L4ROPhE5uVDca9275SIxFO0yRx+TGPkNX9lObRhFrXVgLb9TrKtFqywxE6YmhWO+fvms/NP97MZf0u445hdxzU7x0KDW3Ia4pz/AhQAnxG7WwV5c1YRqANnONdSwm9eQpaoggJ+z7ES8Ch0XBZuyzsWi0zrIPoOen9g5b7Zus3PL/4CYpCLkT8uksh6BgK8aDdwWiff/+TBkyC814/aE0g1qkXTmv68XHNxVvtXPbWUrpmCU4etZmPN39EZbCSHik9ODX3VAZnDibTkok/7Cffkc+8nb8wZ8cPIIJ0MmVzUXkJZ5bujA3MdRAElpuMLLSY+Sk5h23hKEIa6Zk4hJGp5/FFXogSV4C3rziaUd2aEK6RNwN+eKTum26ch9NT+cFi5uedu9HWd20bcY6jMsqgGSPJ0o5gTrvUg7q2++ILRpj85hJ+2+ngryM60zPbxndri/hlSxnH9cxg+oWD926eq8PO3Tot92aks9Jk5ByXm3vtDixRCQJE1gA441+1bhS73buZ+sNUdrp28tQxT1Fp78fdH6/mlH5Z/PviIfvPqtfRhoq1Wt5PtDHPYmanXodRSob5A1zqrGRkj7MPvd3mzYB5T4KnhB8tZu7ITCc7HOaxsnKG+WND0XKTkSfSUtiq13OfsQsXXfTVQUkFIgEu+fpKNjvWMm2Pg7H+mBO175gQAR7MSONrawKPlto5t8vph27nrqVEP72G5b5CNhoMBIQgJxxmlM9PSl19Z8ythzQjF46G+XLuPXy041s21LNSkBUOc5bbwwSXm47hCNTRhgAiUYnDG8QfimA16kgy6+ufDT3IceiQ2KevbNHreTkliXkWM+Ea5dRIyXFeH1M1mfQ5/YWDdqoAQkvf5LHl/+Rzi55jvD7+VlZOu8jeRdcQ8HZSIi+lJNElFOL15JGkT3z7oPWAmCP16bUs8e3mhZRk1piM1V/ppOQkj5ebHRV0iD9UHuq1LfeXM/WLSaz3FXGuy8OllZV0DYUBKNJq+cRm5Z0kGxK4s7yCiS43ovPoQ3JYq+yMfnoNu1y7KNVp0UjICYfJjtSzqN0CbQhidbhdr2e7XkdACIxS0iEcpnswhH6fvhIMRyn3BIlKSZJZT4KxCQv28b6y0mjgibRUNhkN9AsEGO/x0TEcxqURLDWZmJNgwSQl95i7M+HCLw7NzngbwrGVCPCFNYHXkxMp0O/dd9UjGOR6h5MTvb467ysHwtaKrVz81UTaB8P8n7YnmePvOrS2cRAcqnO8rY6PpZSyaxOETwVeILYD6A0p5dMNHd/qzvHXtyHz3mLfobzqilR9XqjTMrldFgDve4xkT3itSZVY5g4wb2MJoVdfYrdnDu+cEqJfIMB1Dicj/AE0SH42m3kpJZmtBj3XO5xcX+HcrzykdIVzXz2ghlMV9G75+koWFC9hlclIqVaLBmgXDtMlFKZ3IEiXUIi6thz6bblM9k5iZ0YZ0rYIX9jLuI7juKzfZQzJHFLnDdBut7N97Tx+WvYUeTYXa0xGDFHJWJ+PQf4A7cJhIkKwQ69jhdHIbyYjPo0GvZQM8QXIDOtx5x7N4sr1hKNh7hr6IG/NTqPcE+S7244j3WqsUxMgbf69sOZ/jV6X83KyyYhE+L/i0th1ruvaNuIc2+12Tv9oCsmGrczeXVf3aARTCpz4SK0n74c+X8t7S3bw0kVDOH1gOwAC4QAf5xXxyFfryDGFmX7RYAaueqyWnVHgU1sC/0xNQQB/KyvnNI+X/RFwxrRams6Ak5t/vJkVJSu45+h7CJaP5tGv1nPu4PY8O/EoHI7Y82/arKth64/V54WBd5JsvJycRFgIRvr89AkG8QgNPySYKdXpmOByc18kCfO5TesrVdjtdtj9G2mLHoWiNQAsMxm5JjuTPoEgrxSXYIvKKosA8AvBXRlpzE+wcK1XMPW019F0HtEkvXAkyn9/Wc/zK58gkLyOZ4rLONXrrXtMiD9Ah4CbsjJYbDbxf0UljDK3P+j+mTrvHr7ZNovpKcns1te+aeqk5HS3hxsdzv1v/tYsOP7+Rmdv3IEw6wsr2VrqpqiklJJgPqsqn2OnCNA7EOQct5tRPj/tw2EiCIp0WlYYjcxNsLDQbEICx/r8/LXSxShfgOCpz7Eo+Ux+2VLG4q12tpS4CYT3OvAJBi0DOyRzTI90zhjYDiuxB/6abWiHTsdPFjN5JiPFOi1BIbBEJZ1DYfoHAoz3+siqsreOvtLUa1tzTAgCL6Qm836iDUtUMsHt5jivj4xIBLtWy0KziU9sVpwaDRNdHkakXk1x17+Sk2ymS3oCnVItDYZhVWnqv76CW3wbyDObuC4+nledte995VeTkZuzMsgOR3jTZyKzifeVWpq7fyNt+b+I7ljEtJRk/pOcSPtQmIsqXQwMBPALwQKLmU9sVqLA/XYHE9zxmc+DuK9s21XEvK2b+Gjr7ZTj558ldsb5fLWOkXEbi7RaHspIZbHZzHiPl0fL7CRF2W8cqj5PSkpdATYVu9hS4iYUiWLWa7EJP72yrLiW3s5nZXn8bDZRqa39UNcuHGacx8d5Ljc9Q6HaP9xMbUgCS01GPrZZWWAx49Hs3x4s0SgjfH7OdHsJZd3Ec2Wj2W73Eonu9bPSrQaO6pDM8C6pjOmWTt+cRLSaGtmeZl1N+fb5PJ+azOc2K9nhMPfYHYz3+vYbl3bpdDycnsoys4mJ7iBT+txC9rjrmxyuU63524vVD65+IbgnI40fEywM9Ac42+2hUyjE7wY9H9usbDEYOMXt4fGycsySeuuzIc1yr5Mbf5iIX/r4sKCIdtEIQqOHK75pVQf5kJzjQxDVAvnASUABsAy4SEpZ7zb/1neOb0Xm/Wd/ZzROzcFsk0HPZe2yaB8K8589xdhG34w46bH9zilweJm9tohv1uxhxa4KpPDx6QdXENJEyLs5nWsqnNXOaNXvB4TgibQUvrBZubDSxX12R91vZ2lkxmhXuZdZa/awcmcFWetfIyq+ZW53Gy6tBmM0SlYkQgRBsU5bPWtiiEq6h4L0DoboEYwNKnviN8g1JiMa4JTcv3D1wKsb3ZiY/8nTMP8JeqRpQUK+Uc9Mm5UFZjOFNW78Qkq6hUIM8wc41utjmD+AWe51eMpGXse9GjtL9izhzkFP8ORMA2cOzOG5SUfVrbnk/+hpqdjvu33xCsGozh2YUlHJjRX7vNem5rVtxDnO/+lj3pp7H990NbJoRwHWg+1D8VmNFTsdTHh5EVcd04VJo7W8tvo1Fu9ZTGWwEpPWRDtLF4p/83BucCMTkpykRCM4NFryTEZmJsYGq6N9fh4vs9MuHKm+O9XZrvdpQ/6wn3sX3MvcnXO5qv9VaCtO47k5m5k8sjOXJG9AzPkbPfWF1cdvMuh5KD2NDUYDJ3i83FXuoEM4Ut2WQ8ArKUm8kZTI4ECA6UWlJIy6Ge3J+/eVOq9tvA31TIvd/Iq1Wia1zyYpEuW9wiJsMrZFTlTdgeP/hIG/pafylc3KdQ4n/XSnI098lFHd0uqMoy5x+flyZSEzFm3HUPEYBcmF3Kf3cnmlq9HxAMAjBJNzsijW6nhvT1Fs1uwAZnTzf/oY93cP8n77IN9bE+gbCHC508VInx+zlGzR6/namsDHtgR0wO1Vs2/7/lAdM2PBcJRv1+7h4+UFLN5qJxSRQJh08QaupI10Sddxh93Bib7YzXbf1lulUTUDONNmxa7T0iUY4qJKFyWOY/k3lzE8N4V+OUm0TzZj0mtw+cPsLPeyfIeDdYWVAJyoWcHVwbch3cnPFjMLzGa2G2IjYG4wRIdwGJOUuDWCrXo9JTodQkqO8fmZUuFkcCBYr50NXVu+f6i63e7SabkzM531RiMTK13c7HCSHI3WslsAlRrBy8nJfJBoJTcUZnhhX970XQWAQaehZ5aV3tmJ9M620addIr2ybdUP7BvnzaR0zkM830WyTa/n8VI7Z3i8yH1Fqog/aC03GpmanUFWOMKbe4pJG3kTmgPsK13StPwtvpoxqdLFXeUVmPYZk4q0Wh7MSGOJ2cTFThd3lTv2bjhqpN1W+kPMWVfMrDV72P7bDFxZ36NvZ2B6cSlDA7EVnJpykr0rLlHgvUQb01KTyYhEeKakjKMCQfx9zmPdiGfJL3azqaiSTcUuNhW5cHhD++l39H2KJfUnCnPN2CJRxnu9DPYHyIlPtmzX61hmMvGzxUxICI73eLnR4aTXvk7yQbYhCSwym3glOYmVJiOpkQjjPD6G+f10DYWwRCU+jWCHXk+eycg8i5kSnY60cISjgp3IyH2ELqnt0AiBwxtkW6mH5TscbC2LPaQkW/SM7pbGUN9Shv8+jQ3pbl5MScKr0XCps5JrKiqx7HuPif+vELHVrH+lJPNOciLHeXycXNaLFYOe4eIRneiZZWuynQBOjYabstJZaTRyT7mDiyvdtcaIMLFVjxdTkugTDPLvolLSo9Emj33RqOT9d1/l423/YmeugTf2lFS3IQDG/w2Obb3wikOdObYAtwOdpJTXCCF6EHud9NeNnDcKeERKeUr8/+8DkFI+Vd85bRFWwYzTIRIfgA0J0O6o2NJP8fpaMU0S+NVs4oasDIb7/LxUXEq+eRhf9J+OEIJSV4DVBRX8Xhpr8H3bJTK6l4YlpXfy0MMr6RAOk3Pl3k1DVVfdFTVTnHY0nc66l5eWP8sMz2YmVrp4yO6o+yZdx1LGTruXf8zeyLdr9xCV8I/Ef/NOeD2bjQZOsUW4pNLFEH+gejCMmFPY2XEoG3qMZWPxCjbuWcpGraQi/jRuiEr6BQOM8fo52+0hu/NYuPSzhq/lJ1PIn/8hQLVjU9NOu96EPac/2pE3krVzGQl5MxAhT72OXDClC1d17sJGdwHH2x7nk8URZt96XO3OXo9mNVVxmt1Pgi3fs3z3Ei5PtzC9qITj6wphqbq2k++O/X9dzvGch8n/8jnyjEaeGpDDG3uKGeGv0blrxrttnFUrzrEuZNcTuNh/D/nFLm49t4znlv+DBH0CJ3U+iXYJ7agIVLDh99nkbd9FUKPBmF179rxHMMgUd5BTskfxTeIkdi6ayaWGeSRIT71O3r5tKBKN8OSSJ5mZP5OTO5+M1T2JrF9fZHxlbJmuV5oWnxC8mpzIjKREkqJRHnC4OSm1L2UjHmDpkgUM3/4KGTirK/K7BAv3ZaQx0B/g1eISHOljyJ46q+EZjfi1hVh9hoAr22WxyaDnw8Ii2oU0rI505Yf213N56jpyfv8vmhrXNkIsJvjzuIM81p7Mk+JqjF1GkpuWgFGvodIXYn1hJWt2O4lKuDznWT7w7eYMl4fndHU4nzVjfeNjgoyPCYVaLX/NycYso7xXWExaNFpv+EE4EmVtYSVLttrpv+Qu2D2XJ9NTKOho5TZHBZc59zrlEY0B2g9Fd/Jj7Fr7Px7d+RVLTAaO8/p4tNQeuxnVJD4DWJp8FB8s2cl7S3ZQ6grQKdXCqf2z6Zrj5eO1t7B6dwknebw8o3WTUBXWFf9PpTTza7QfizSDGa9bTb/oJtKoABmbdf3eauH9JBvrjEas0ShnR0ycPuouevc5r850l7/t3sLP39/Lqu1L+M1kJJpjxhCVDPf7Gev1cZzPR3sMe/tK/NpuD9iZZU3gv4lWHFotI3x+bi2voF8wyC8cxY3iQVITDHRMtdAn28bw3FSGdk4hJcFQZxuaazHzUHosJOuJMjvHh3XYM0bynu4cErbOZqKYS7Jm70qLABabjNyXkU6lRsNNfiMduz3CknB3NuypZGORi1LX3v5uMWi5TXxAV+eXPJyeRqi9iWnFpYysMSZEtGZESke0I+N7HGq0IYA8k5GpWRl0DoV5o6gEn7EbFeOeptew8fVvyI3b6UXwSt8sFlrM3OioYEpFJQIIomOd6Mm7CVdwkiaPce5Z6KJu/pWazHtJiRzr9fHPkrLqdrDvLHKFN8iPG0v4Zk0RP+eXEoxEuSblLWZGVmGUkndw1pqhlUBUb6M0YwQ/pl2Ir2ANZzpmkCFj/XOt0cDdmWns0em4weHkyopKtkezuCt8PZv0feiZbaNXlo1e8X97ZNnYWrmOV+bfwcKC3aSFI9xm8nOW24NJ7s0hE9UY0HaM7XNw7M7jf8un87ZZg1sjODXuJHcKh6vLGU7uyprhT7Mw0JWVu5ysK3TiCYQRQpBmNdA+2W2IsC4AACAASURBVMy1oXfIWPc2EijukMArKUmsNRrJDoe5qqKSCQGBMT4m+FZ/iWbFDAxhV/XNLiJgodnEfxNtLLCY0Uk4JXsklwy7hf7p/avLUuLy8+vvdn7ZXEb/Dc+hd87i7eRECjtZGe7zc7+9nG6hcHX/LJapfGw4k44pFiY43yYhZK/Vbj+yWfl7WgqDAwEu2p3JdcH7Gdk1lctH53Jin6w6Q+Vq9pUirZbrsjPYqdfzdEkZJwYilCYO4FXDpXQo/pEJ0Tkki1hf+cli4u7MdDIjEV4tKiEnHnol6gmzCIajzF5XhG/Wg6zwzObTRCtP6txMdNWI3xYauPK7I2fmWAjxX2A5cKmUsr8Qwgz8KqUc1Mh55wOnSimvjv//ZGCElPLGfY67BrgGoFOnTkN37NjRRLOaiV1LYfsCyD227krZtRQ+uxbKtwLwmTWBv2Wkcbrbw+OldkJRA+/JU/iP5Qr65SQyPDeVU/tn41n/EjdsfpcAks8fzscajWKZUjt1mG/4TdzvPp/PVuymV5aNZycOZG7xDN5c+yYXBgT3F+6o38EZcyuhEx7m9QVbeeGHzeg0gnsHuOhV9Dh3W8KUlwS5o7yCyxJC1QNJedTKG4ZLyBx3HZOGdawV++T/fRG7Pp5Chr+QBBndP9SivqWpGvFY+fbYUmiVoyqJzSh8pz0Ow8Q3GNcrEyFg2XYHj361jgFFn/GQ+X9YInXP1pVpNJyXm0t2Sm/W5V3KcT0zefmvQ2vFRu2rWfP67PskO2PtDJ5b/hwzCyP0Duyu78rCZymQklvbOa6xSSHfHsElBJcPzeWm8gqucVY2vHS3TxuqiQS8UQNvdRnHa5oNjO0wliePeTK2c7eG5jp7hF16PdF2JlxCg1VG6RMM0vXoGykecR9PfbOBz1cWcnLfLKZfPBjjnuUw6/bq0IQ6qXGNpJT8Z91/mP7bi5gjEc6trCRjlwchwN7Bwuc2K+VaLRO8QW4ffAtbcy7i7UXbmbV6DxI4fUA7buxRTo+FdyIcMTu/TbBwd2Y6p7g9/KPETkgYcA+6ivRz9omu2ufaQqw+n05N4f0kG8+UlLHNeQIfJl7FP84dyOju6fVe2yixuPLPbVaujS9rz9Mfzy2B6wlFJQkGLb2ybZydVoBhz+M8YpMM2VrBveUO+jShDdXUlOVbWW00cFV2Jn2DQV7fU4IhXqfbel7N6j63scPuZcXOCvK2l3Nm+Hvu1H7EHnOIKTIRr0YwXePiuPiDmgS+ih7DLcGpmPQajuuRwan9szmhVybf/HIXz+2ZizUa5bFSO2NrnFP17xfhMdwWvoHje2Vw+ehcjuuRwazlL/LEujfQRyNM2VDCKL9/v74SHn0Lm/rfwapdTraUuKnwBvEGI/SLbmRi8fNkeTdXH7vaaOD9RBtzEiyEhcAstHRJ7UWaKQ0hBBX+CnZU/I4zHLvpWQu8DPMHONcW5mifH0tUUiGsfJ0+hdKeF5FmNSKlxOkLs93uQVe4jOvLnyVLW8wniVbeSE7EodVyanV4iZb5SROYrr2ETUWu+Kw43Jy0kKvCH5IYKSffHiEEfNkjnQ+TbPQPBHimpIzfky/gbud5lLgCJJp0nHFUDmcdlcNw7Ra0395Rq6/YNRoezEjjF4uZMV4fj3c5l4xTnwFi4XKbilzYNy5g9OoHydOXc084AYOUzMBJ32AoNp4dwJiwyGzixqwMegWDvFZYQoKUrBB9mddxKpZuo2ifbMas15KR/yE910/HErKztFzyWEYqBR2tPFRWznkuD1LAUtuJfNHlYYJhidMXYofdw7YyDwPlJv6pfYW8ZDdPpafQIxjipeJSMuMhLBJYk3oKTxpvY/kOB5GoJDvRxDVdyuhY+gT326Ik7PbxeJmdY5P3saeOvhKNSgrX/kTCnLtJdm3CJQSPp6cy25rAUJ+fmxxOhgQCyNG3VM+WSylZXrycGXnT+Mm+itRIhNM32znN7aF/mrZWe/8sPIY7wjcwqGMyJ/XNYnCnZHpl2dAULWTG3Fv5wBAhLATnudxcV+EkLRIFGT83MoaXU+5mYIdkksx6pJTYSldwTuFz5Ea28a5Hz/+SbBR0spITCnO108lZLg/5uVewqMtN7K7wsWa3k1W7KohKODVpB49o3yLLu7nWvWyHTsdHiVY+s1nxaDT0TevLMe2PoXdqb5KNybgKf2P9sn8zWxMkvzxKVijC/UYvp3hioV0eTSKLu9yIb+BkhnROISfJtHeCoY77yrcJFu7PSKN7MMQLpW7yI3150XcaxUlHMXlUZy4c3pHk9R/sN96KbBPXZWfg1WiYVlxKgDFc7ryaqITcNAsju6bRLcNKrm8tg9c8QZp7E6uMBqZmZZIgo7y6p4Qu4TAS+D37dDaPeQ6NEJS6/Py2swLHpl94KDidtcku7opYOc3tYZrWXeNK7R/21xocqnOcJ6UcJoRYIaUcHP9slZRy/zXu2udNBE7Zxzk+Wkp5U33nHNavj/5kSnXs0RtJibyQmsyxXh9Pltpjm2Y0OjDakAi+0Ef4e4qN5GiUl4tKyXk5Fshf7RzvM7v048Zi7v1kDeWeIDeO604w6Qve3fAOl+izuTt/aZ2OowR8GPFFdFi0UXRawacWDU+npdA+HObWdUV0Coerb4TRricwd9grvPbz7yzb7sBq1DGiSyo5yWbKPUF+3lyKyx9mZsZ/GOaaU79TrjWCVg8IiAT2zrrD/o5q1gA2DnuEqT/p2FrmITXBgEZAmTtIVqKRx8/uz8n9smtd2335OsHCfZnp3GN3c1aFB4tei65GBoG6NOvbJHDjD7cxb1sed/Z9hyu8M+rfHDTDE3uKnZLRoJ13DuxAp1CIl6wDG59Zhzo3JElgrcHApTlZDPMHebncix4J0SiE97ezR9zOrZEs7pc3UGAdQKHTh0YIbhnfg6nHd6s9Q9DYJqiq+tTqIRImX/p5JSWJuRYzvuKYzQlZBkb5/JwSaM/POS+yYmdsSdBq1DFpWEeuPCaXDik1UqnVqM8ZiTaeS0vhUmcld5bHZiJDQo9Ob4ilIarn2m7pZOOezHTOcAo27bmC3kefyAOn9al/M0sNzSixGeTPbFamVDi50eFEo9HFZvXj1/ZDc6yvHO33c8u6Igw0rQ3VpTk7wcJdmemc7vbwVKm92mP1SR1htAiNFiNB9ISZnWDhofRUzLv9PGQv5+Sk+CxwfOYu2G4YedvL+W5dEbPXFVFcGUCnEYzqlkbvTl4WF97PdtxMqnRxR3lF9XKrjK/ARIUOrd5MqQaeSbIw22JkiN/PP0rsVJbErnO1nU2NOa2jDZVrNCwzGfnNZGKHwUC5VosQAlskSqegnx7BECP8fkJFAURcUwKlWWN4KvXvrNjpYEd57bCDnCQTXTOs9M62cV75G/Te+iZuIXg7KZF3kmyEhOB8l5vrHE7S0BA12AhHo4igFx2xWUwBfFWp5aWUJAo7WZlcUclZZUbuD1/PZmNfxnRL55zBOYzrnVmdlaU+OyXwX5uV51KT0QGTXV7O8UXIiUSIRqOs1ob5T1IiPyZYyN3h4u5yx16nsesJTRsTarTb+WYzt2WlMzAQ4OWiUizxGVJf1EAYLQaCmESsfxTodUwmiTKthhe1bsZ5fcisAWgamLnbUuJm7W4nPRbejjvwM3dmpmONRnm5uJRewVB1XYSElojGiF4DGuB/8b7SMxjijrV7SJHRg25DEvjcmsC01GTKtVq6BYMM9gdIlBrK9DqWG3Ts1mmxRaJc4azkr5UuCspiM7/7ahZY+/P16j18ubKQ9Xsq95O8wfA2/vSlfJpoRSclk1xuzne5yY1vHBQ1xwQEjrCH2QmxFZJ8e5R2oTC3mfyc5vawO5rNHaHr+E3GQgutRh29s22M7pbGiX2zGNA+Kea01jPeeoTgc1sC31htrDXoiO6zGXRQIMCIrU6O83rpX2VnU9vQPpq/mE3cnplORo1Z3QAGgtF4G9LEHWJi4+0ag4Hn+rXDJKM8UBjl397rqEgfzOkD2vGX/u3o0862/4pfXHOjQc+12ZkAvFJUQu9gKLbaJLX4MSCIotVILAT51WTihuwMuu1w8WiZnX4HOt62AIfqHC8CxgMLpZRDhBDdgA+llA1ackSEVRwoNWZJ/2ez8lRaConRKJMq3QwIBCiLx+itNhkZ6vPzTKm9+qm8mnpmoyq8QR7+ch1frCykf/tEevSeyw+Fn3B52jBuz/t0/w1CNWIuQ8CT6al8arNyrNfH06VlJNbYALCv5m87HczMK2DZ9nJKXQGSLXqGdErhryM6MSw3tUlZHxqlhmYgHGHW6j0s3VZOJCoZ0jmFswflYDHUcHTqmV2VwOXtMtml0zFr1x5MyAZn0xuKexr74UkUl2XwzhkvcXSXuJ1f37L/gTPiSz2XJzRo4oPpqfxgSeTXK1Y2PV9pDTslEBRwQU42Lo2Gz3bvqV1v9TFgEr8Nf4ZvVu/B7gmSm5bAhMHt6ZRWT67f+uxsAI8QbDHoERK6hEIs04/l9vANmPWxmdcT+2RxzuD2WOtzVuNtSPodPJ2awgdJNu4rK+cil3vvjp04+165jQY9k9tlkRKwUlH+GP88fwhje2Y0XugadkaBR+N9YqjPz9XOSvoFguzU63gjKZH5CRbGen08U1JWO57vQDNBxGe9X9P5mJ6azFkuN4+UlddKIC+I9dGqJe3Bfj/PF5fFwjAa0IxGJasKKpi9rojv1hax3e4FEaZjxkdUpK2lfSjMFc5KTvZ4SYnH0G7T6/jSmsBHiTZCCKY4nVxdUbl/QvsD3cV/EG1oP/bRrMp0oRECi0G7f2x4jb5SqtXwanISH9usGKXkUqeL0zweckPh6pjzVUYj/0208q01gfRwhL+V2cm2jmP9yGfpnmmlX05i47nN67Bzu07H86nJ/JgQ61+2SJSggIBGQ0I0ylUVlVzurNy72nYw1zY+kzc7wcI9GWkM8Qd4rqSM1Dpio5eZjNyWmY4EpheXMSQQOPB2u2sp6z+/hpssAdwaDQ+WlXNGfCNvVX90ajT8MzWZL2xWxnp9PF1SVnt/xYHaWaM+fULwhTWBHxIs5Bv0uDQaUiIR+gWCnOT1caLHu3+cbQOaDk+QVQUVbCvz4PCG0GsEiWY9g8jHuuwuXtO5+D6+2pEbDDE4ECA7HEEgKddq2WAwsNZoICIE/eLx/yd6vLF+M2AS8tzX8AQjeANhLEYdCQZt/WN+I33FLQS79TocGg3WqKRLKLQ3vKURO+tln/vnyvisrllGebWolO6hvSvINUv9qTWBx9NT6RAKc7mjO0XDpzG+Tya9supwiOuxc4dOx5R2mbiFhn8XlzAoENzv0F/NJm7LTKdDKMyMPcXYquxtjmwih8ChOscnAw8AfYHvgTHAFVLKeY2cpyO2IW88sJvYhryLpZT7J8iNc9g7x1W8MwG2/sgmvZ7nU5NZZDYh4w2pfSjMlc5KznO5q3MXu5cGAIH1sf9rdNngmzV7+NsX6yhz+2nfbTaVhp8Yn3Ac929bQoZnc61AfICdOh13Zaax3mhkSoWTGxxOtIDdG4XkXNImv3nwT2RxO5tKs2jWMYucZzJyRbss7rQ7uKzSdVCaZb4yxv1vHIGS01hx85N7ZyB3Ld0//KAR59jujUJCJtNzj+cTlvLZmV/RPTX3wOx8ZwJy64+8npTI9NRkXi4q4di64qBraiJIm/TiwS091WVnI9i9UTAmk3b24we/3PXOBCJbf+S2zHR+sph5vqSME7y+Og+VwO9+mJKdhduayKjEZ3n8jDF7U9g1harE9TsWVc9S/Ss1uTqeHmI7yq+rcHKp09VsfUV+fDWvFHzHyynJdA8GucnhZLTPTwT4yWLm/5KT2G7Qc4mzktvLK6g8CE2nL0RhhQ8p4Xf3St5bcjfro/FNPZEIQSHwajRopOREj5db9om3PJLbUNU4tF2n48XUZObEHdXESARrVFKu1eCPO6un74nVbeeJLzSrnbt0WuZbLOzU6zBFJb2CQcZ6fdikbB4742PfNwkWHkpPIyka4QaHk1M8XqxSsl2n450kGx/brHQJhXlsWzEdwpGDr0+geOZl3F7+K6tNRnoEg4z3+EiORtii1zPbmoBPCK50VjbvfaWBlcK6aK5rW7buY76xJrDYbGK9wUC5VoMUAms0Ss9gkKH+ACd7vPQKhij/A/SVfL2e67IzcGs0THU4Od/lrn64qepH30oDw30Bph1zJ6kjpx64Xny8LSpYwpTsTAp1OqY4nVxU6SYpGsUtBP9JSuSt5ES6BkP8X3EpWlfo0O1sJg45W4UQIg0YSXy/gpSyrInCpwHTiKVye0tK+WRDxx8xzjHUml0t02go0OtIjEZjyfP3OdT7joDE9lhW1vtcUAtfMMKHS3cyc/lOtsp3MaQsIezpymmubB4LzMISdVGpEcy02Xg9ORGdlDxWVs74Gk5HftopcMoTtZN7H6KdjZHfeTKMufnQNeuYRZ6SncEWvYHZO3djqLGBr6l2ztkxh9vn305a5e3Mv+mK/Q+oaWcjznGVnQscJUzbeD1T+z3E9cMmHZCJ0ajksadv49t2PzDC5+PFkoa7VL7xKBh9Mz3Hnn9AOvtxIPXZjG3I98MjXJViIN+gZ1pJGWNqPgjEhyCnRjBBk02RQfDchBmc1rNp6djqpEYbCgKLzGYK9DrSIxFG+3y1Zuib0875Cx7nKZuhVnYWiGVmuKvcUR1f3Bx9RUrJ+vUz+fWXJ9kTdGKQki6hMMd5fXXmfT3S21BNzV06Lb+azeQb9Hjis44DA0GO8fooNAxsPjvre9HEPjSrnfP+zsZQOY+lpcYyBkmJRUrcGg1aKbmw0s1Njgp2N5Od0Z2L+XLWtczU+Fgdz49sjkY5wevjiorKWlkf2uS+0lya+9xXQsTuI3WtfeXTHcbdd+T2lXi73aPV8kR6Kj9bzJji/olPCLYb9BijUc4uMnLeyHvoO+6CQ9PctZTyz6/hKU0ls60J6KQkKxyhRKclJARnuD08UFaOVcrms7MZONSZ47lSyvGNfdYcHFHOcRV5M2DBc+ArrxGLG4xlvugwHMbcgvf8WJi1ZcmSA/75Sl+QjzZ+yn82vIg75CLTkklC0MeuUCVhITjeH+b+CjftIjL2KuisvnDio+T7YsFvzdYA5zwMK96DcIC96+Kylp3Nrlnj2i4y6rk23caDpZWc4QmT0GngAdn51JKn+GD9TE5KeJ1/TWrghY9zHobrngEZhavTG7RTn9qOM78cx5C08bxz9j8OzLTt5fz10/sxpf3KZ90vpevSGbE2VKWn1YMxEbIHtMy13bc+q9qujLZYGypf/DLXrX+FzZootzm9XOIJokES0SewIbMf95u9bN6+kwdGPMAlYw/ubZT7UbN/ViclEi3aV0LfP8Ti9R+xVivRScmgYIQhwRDaVuorrWVnW7Shpoy3f4hrmzcDufjfLA+Ws0QvqNQIOoUjnBzRkZHTcnYGFjyHK+AgNQqatmhDrXxfqTXeao3Vr4POTxzdvJpVdsoI1W22FdoQi19GekpZrYkw22Jkm06DRWjpZ2nPWSPvwhHNbXbN9YueY47Gxx6thsxIlFN8IfpFaDk7D4GDco6FECbAAswDjmfvZF0i8K2Usk9zF/SIdI6bQM3XRx8s7qCbr7d+zcrSlQQjQTrZOnFS7kn0S+tX5/F1vhayhWlJTSkl5355LqWuECUbp7LsgZNIMuubrHnO5+ezqTDE/YNfYPKo3IbFGstzHNfs0aMHA1+ZhNXq5NfJ3x+ANXDPZ4uYVXEjp3c9lX+M/Xujx/9R6rMyWMkDvzzA/F3z6ZLUheFZw7H77fxU8BM2vY2p7acyJGvIEW+n0lSaSlNpKs3DQ7M+GnKOG3qP4bXArUAOsVRuVc5xJfDvZi2holGsBisX9r6QC3tf2NZFaROEEFza91L+tuhvRAyb+WZNfy46ulOTznUFXWx15hPxnsDADvvmIDq0MmXp+1Mc/ZhiTzFZCVlNOi8UiTJ718eIpBDXHnV1s5XnSCDRkMiL417ku+3f8d9N/2X29tnYDDYu6HUBV/W/CkfBIWwCVSgUCoWiGWhKWMVNUsrprVGYP+rMsaJ5CEQCnDTzJDyVHejBTXx0zagmnbegYAFT504luGsKq++9Yf8UTvvSyMxxTW77bBY/VN7Lk2P+zlndz2xSeWav386diycyIG0oH579SpPOUSgUCoVC0Xw0NHPcSF4bkFJOF0KMFkJcLIS4tOqv+Yv5x8X90ku4X3qpVTXtdnv1e9P/KJpGrZHze55PyLiOpQVb2OP0NUlzefFykFq6J/Vv3DFuAjU1x3QagIyYmbt9YZPPf23F+witnztHXHdQmq2F0lSaSlNpKk2leaRrHgyNOsdCiHeBZ4FjgOHxvzo9bUXdaN59F82777aq5h+10U/qNQkhBPrkxXy9ak+TNJcWLUUGOjC4YxNy5TaBmpoDO6QQ9nRnafEiojLayJng9PnI939DmrYvQ7MbfMlkvZqthdJUmkpTaSpNpXmkax4MDcUcVzEM6CubkvNNoWhhshOyGd/pBOZGF/LZqm2Mze7Y4PHl/nLWlq0l6BrP8BGpzV6e7hlWDMEBuMNrWFu2loEZAxs8/l+/foDQVXJJ7zrSySkUCoVCoWhzGp05BtYC2S1dEIWiqVzc+2Kiwstmz8/sLPc0eOzC3QuRSMLu3ozsmtbsZdFoBEPTR4HUMG9Xg+/FIRKN8M2uDxHBDlw++ORmL4tCoVAoFIpDpynOcTqwXgjxnRDiy6q/li6YQlEfQ7OG0i2pB/qURcxd33CC/gUFC9DLRDpZe5CVaGqR8hzTrTNhby7fb5/b4HFfbJ6Nn2JGpk5s/DW2CoVCoVAo2oSmZKsYW9fnUsqfmrswKluFoql8uvlTHl70MEkVN7Pg5qvrfA98KBJi7P/G4rL35vR2t/L0eQ2HPFRzANkqANbudjLhvX9iyv6C/53xP/qk7Z8CXErJSf+dwJ5KJx+e9gkDOzR/iIdCoVAoFIqmcajZKn6q+QeEgQN7V+6fHJWtovk5rctpmLU2in3fsmDttjqPWbB7Aa6gC6+jf7OGVOxrZ592iZiDw9Cg59PNn9Z5zvc7vqc48DspoVMZ0D7lkDVbA6WpNJWm0lSaSvNI1zwYmrS2K4QYJIR4RgixHXgC2NCipfqDobJVND8mnYlzu58HrOPtBXUvYnz1+1eYNIlIb3fG9myeTBWwv51ajeDEnrlE3QP54vcvcPhrv8giGAny7LLnifizuajvhDpnuQ9UszVQmkpTaSpNpak0j3TNg6Fe51gI0VMI8TchxAbgJWAXsTCMcVLK1p0GVSjqYMpRl6MVBhaWfEogHKn1XaG7kHm75qHzHs2wzhmkJBhatCwn98vGU3wcvrCPN9a8Ueu7V1a9QpF3N8GS05gwuOHsGgqFQqFQKNqWhmaONwLjgTOllMfE35IXaeB4haJVSTOncWz2X4iaN/HSojm1vntr7VtIoGjXME7u17TXOh8KY3tmYJDt6GwYx3sb3uPXwl8B+Hbbt7yx5g2Eezjjc48lJ9nc4mVRKBQKhUJx8DTkHJ8HFAHzhBCvCyHGAwe+HqxQtCA3DL0YbSSFd7c8jd0XW6pZXbqamfkz6WYajyaawtmD2rd4OcwGLSf0zqRgy4nkJnZh6typXPD1Bdz98920N/elsuBMrhzTpcXLoVAoFAqF4tBoSraKBOAc4CLgBOBt4DMp5ffNXRiVrUJxMPzrpzm8tfVuMs3ZTOx9Fh9s/ACT1kzxxqmM7NyB1y49wBc6HmC2iioW/V7Gxa8v4eGzcynRfs0mxyYGZwzlw++7kZZg5csbxxxUvLFCoVAoFIrm5VCzVXiklO9LKc8AOgArgXubuYx/aFS2ipbV/GufQaS6plLm8fPyqpfpYO3ACPN9VHr03HRCjxbRrMvOUV3TGNghif/7cQ9TB97OW6e8ha9kPIUVEe77S+9Dcoz/TPWpNJWm0lSaSlNptiUH9CYCKWW5lPJVKeUJLVWgPyIqW0XLajorHEw753x8W+8k1z2No42P8P5CLxOHdmBAh6QW0azLTiEEj53dnzJ3gKtm5PHUNxt4ef7vTBrWgdHd01tEsyVRmkpTaSpNpak0j3TNg0G9pkvxh2BIpxSmXzSUHaVhXpi7hbE9M3js7P6tXo5BHZN5/oJBrN5dwas/b2XC4PZtUg6FQqFQKBQHh66tC6BQNBen9s9mbM8MHN5gm2aFOHtQe07onYk3GGmxV1YrFAqFQqFoGZRzrPhDYTZoMRvaPl2azaTHZtK3dTEUCoVCoVAcII1mq2hNhBClwI42kE4HytpAV9Ewql4OP1SdHJ6oejn8UHVyeKLq5fCjreqks5SyztfnHlbOcVshhMirL52Hou1Q9XL4oerk8ETVy+GHqpPDE1Uvhx+HY52oDXkKhUKhUCgUCkUc5RwrFAqFQqFQKBRxlHMc47W2LoCiTlS9HH6oOjk8UfVy+KHq5PBE1cvhx2FXJyrmWKFQKBQKhUKhiKNmjhUKhUKhUCgUijjKOVYoFAqFQqFQKOIo5/j/2Tvv+KqK9A8/55zbb3qnBQgldBFQuoKK+rOsDRsWLAgqWFdF11XUde3ormJvKBZQQVRcFVAQpfcOAYEQIL3c3u/8/kghCUkIpBB0ns8n63LvOff7zjkzc94z8847EolEIpFIJBJJGdI5lkgkEolEIpFIypDOsUQikUgkEolEUoZ0jiUSiUQikUgkkjKkcyyRSCQSiUQikZSha6ofVhTlA+AiIE8I0as+5yQkJIgOHTo0lUkNYnv+dgC6J3ZvNk2v1wuAyWSSmlJTakpNqdmI7CzcCUB6fHqzaVZHakpNqXniWLt2bYEQIrGm75rMOQamA9OAj+t7QocOHVizZk2TGdQQBr45EICVd6xsNs2MjAwAunbtKjWlptSUmlKzERkxfQQAi29a3Gya1ZGaUlNqnjgURcms7bsm/QszowAAIABJREFUC6sQQiwBiprq9yUSiUQiOR7mrj/IzhwHuXYvcpdYiURSnRMec6woynhFUdYoirImPz//RJsjkUgkkj8xP27J5t5ZG3D6guwtcPHh0n0n2iSJRNLCUJryrVlRlA7AvPrGHA8YMEC01LAKiUQikZzceAMhhj3/C61jzJRYHyUj10mC5xkWPTCCxEjjiTZPIpE0I4qirBVCDKjpuxM+cnyyMO33aUz7fVqzahYWFlJYWCg1pabUlJpSsxH4ck0WBU4/945qS6G3gISoEC5/kE9WZP6pyik1peafTbO5kc5xPZmxeQYzNs9oVs2/SqWXmlJTakrN5tCcuTqL7m0FL2+bxD7bPvbad9G56zK+WJNFXn7Bn6acUlNq/tk0m5smc44VRfkcWA6kK4pyQFGUW5tKSyKRSCSSuvgj38nWQ3YMyV+T584jPTadBHMCOeq35AW2sulAyYk2USKRtBCaMlvFtUKIVkIIvRCirRDi/abSkkgkEomkLn7ckoNm3sMe9yom9JlAlDGK9lHtSbakYEqaz9LdBSfaRIlE0kI4LudYUZR3GtsQiUQikUiait93FRDfehUxxhiu634dAKqiclvvcajmTH7L3HCCLZRIJC2FWrNVKIoSV9s5wEYhRNvGNkZmq5BIJBJJY+Pxhzjl6dmYOj3Nzb1u5r7+91VsAvK/6/7HsM/PxFXUm5/HTqNtrOXEGiuRSJqF481WkQ+sAdZW+ltT9pfU2Ea2dGS2CqkpNaWm1Dw5NVftK0JYtiAQXJh2YZXvLHoLQ1udhaKu48c12xtNsz78Ga6t1JSaf0bqco73ACOEEB0r/aUJIToCuc1kX4tBZquQmlJTakrNk1Nz6e4CDFFbaBeZSpeYLkd8f1nX8wl73PywZVGjadaHP8O1lZpS889IXc7xf4DYWr57oQlskUgkEomk0VmdeQDV8gej2p+DoihHfD+4zSAUdOy0ybhjiURSh3MshHhdCLGxlu9eazqTJBKJRCJpHEJhQYZjLShhzko9q8ZjzDozrY3peHU7yS7xNLOFEomkpSE3AZFIJBLJn5a9BS6ChgyMqpVe8b1qPe60Vqej6m0s2L2lGa2TSCQtkVqzVZwIZLYKiUQikTQm32w4yD9Wj+H0Nj358II3Kj4vz1ax+KbFAOws3MPoeZcwKGo871521wmwVCKRNCd1ZavQNbcxJyvlmSomDZvUbJrlAe/x8fFSU2pKTakpNY+DlVm7UQ1FnNV+SJ3HxYsownYzO8LNF3d8sl9bqSk1/6zUyzlWFCUW6AKYyj8TQixpKqNaIuWZKqRzLDWlptSUmieP5rq81aCDIW0G13lcUVERMb4OlFh2EA6HUdWmjzo82a+t1JSaf1aO2voVRRkHLAF+Ap4s++8TTWuWRCKRSCQNQwjBQc82DEokadFpRz0+PaYH6OysOrCrGayTSCQtlfq8Gt8DnAZkCiFGAqdSukGIRCKRSCQtlgPFHkKGTNpZutWYwq06Q9qeCsAPu5c1tWkSiaQFUx/n2CuE8AIoimIUQuwA0pvWLIlEIpFIGsaa/QfRjHmcmty3XscPbNcFETKzMb/GLKYSieQvwlGzVSiK8jVwM3AvcBZQDOiFEBc0tjEyW4VEIpFIGov7v5vFgqKnef2stzmjXdUFedWzVZRz+vtXE1btrLn5h2ayUiKRnAgalK1CCHFZ2f99QlGURUA08GMj2ndSILNVSE2pKTWl5smlubVgMygK/ZL71FuznaUbGb6vcfldWA3WBunXV/NkvLZSU2r+pRfkASiKoimK0hrYC2wAUprUqhbIjM0zKjJWNBd/lT3TpabUlJpSsyk0c3wZRKhtiDBE1FvzlKQ+oAgW7V3XYP36ajYnUlNqnoyazU19slXcBeQCC4Dvy/7mNbFdEolEIpEcN7k2DyHDPjpE9Dim887q0B+AJfvXNoVZEonkJKA+eY7vAdKFEH/u1wSJRCKR/GlYvGcbiuahf8opx3TegNR2hP3xbC3c3ESWSSSSlk59wiqyAFtTGyKRSCQSSWPxe1bpyO+oTqcf03kmvYZVpJHt3cnRFqxLJJI/J/XJVvE+panbvgd85Z8LIV5ubGNktgqJRCKRNAajPrqb3PByNty0ElU5chyotmwVAKM/e4GdgRksGL2AFOtfbomNRPKXoEHZKoD9ZX+Gsr+/JDJbhdSUmlLzr6AZDIWZu+EQP2/PpcQdINkU4NK+bRjR5+g7zDUWjVHOvEAGMYZONTrGR9M8NakvOw/O4Lf9a7my+4XHbcOxaDYGO3McfL5qPxm5DqxGHWd0TeTK/m0x6bUm06wPUlNqnmzUJ5Xbk81hSEunPFOFdI6lptQ8uTSFECw5sIR5e+ZxwHGABHMCozqM4sKOF6Kp2hHHN4ZmU9HUmlsP2bh35gZ25TlpE2OmVbSJeWuX8fXG/Zw+oB03D+zH0DYDiTZGN4l+OQ0tZ7bdRkh3iE5RQ49L88wOffg8S8fvWetOCuc4GArz/I87eP/3veg1lZ6to8i2eVmwLZe3f/2DqVeewsC0+EbVPBZamubOop28v/l91uatBQF9k/pyS+9b6Bnfs8k0m4q/imZzc1TnWFGUrsADQIfKxwshzmo6sySSUtwBN9sKt+EP+ekc25kkS9KJNklyEmHz2Xhs6WMsylpEgjmB9Nh0dpfsZvHvi/l8++e8cMYLtItqd6LNbBH8siOXOz9dR4zZwNs39Kd1Uj5T105FtS/DCOzwG5n826cYVCOXdbmUiX0nEmuKPdFm18hPu1ajKILTW9dvZ7zq9GkbT9jbmm1FWxrZssbHGwhx28dr+G1XAWMGpvLQeenEWEoneZfuLuCfc7dw3Xsree6KPozu3/YEW3tiCYswb258k3c2vUOEPoIz2p6BqqgsObCEhfsXclvv25jYd2K9thqX/LmpT1jFl8BbwHtAqGnNOUlZMAXWfwIiBCggwpDcA855Etod22KQepG1ipLfX2JOyXZ+1YXYpykYUOilj+ZvPccy4tTbmqZxr5kOv00FTxGggKo1WTkLPAVMWz+Nebu+xke44vN+gTB3WjozcPijTXNtoVnLWUF5HQr6AAEGK7Q9DYbe0zzl1PRgjIKU3s2n2cTXNt+dz7j549hv28sDDj9jDu5CTwbCYOWHVp15xraX63+4njfPeZMe8ceW7qs2NuZv5NPtn7L94AqM7mIGezzc6PSRiNqi69DS3QXc/sk60pMjeW9sf77L/JTJP0wjzhTHuF7jGNhqIEVqGx7+bgFRrTbxVcZsFmQu4JGEIZy/ad7hOqQZwRwDA++AATfVqBUKC9bsK2J/kZuuyZH0aRt9bP1VPfrb5QfXA3B+54H1/91KRJn0WOhInm85wXAQ3bpPYMUb4MqHkL/B7XNHjp0tB+34C4s5pV0to/D1aCuBUJjbPv2FFXnLOWOwE0eEhydX6mgT0Yb+yf0Z2nEocycOZeKn63jwq9ItsftE1mHYie6HmrBPEELw+NLH+eaPb/ibOZWHMncQvXcGILAbrTyflMLbm94m35PPlMFT6h2OU2/WTK9ahzQ9RCTX2VYawn77fraufgPvju9J8Dro4w8RrTRDP1S9nM1Rh5qA+izIWyuE6N8cxrTkBXkD3yztZFfesfLwhwumwIq3IOSt/cTk3nDRy8dVKTIyMgDo2rVr6QdZq2DhFH7OX8+UhDhsmkYvn490fwCPorDGZCRPp6OvMPDvwU+Qmn5xwzWhtLIvfAKXr4TfLGa2GQyEFejm83Om20OkEI1azlXZq3hg4R04Qj4uczgZ4fZgFoLNRgMzoyLJ1um42u5gckQ39Oc81Sia/pCfFUueZuvWmbjCPhKDIfp5ffTy+6ny2G7M+wknpg4t/Q/sWcJu4eV3s4ltRgPZOh0KgpRgiFN8Ps7R4kge+vfj6rRrLeeqdyHgqv3ERiynw+/ghi/P45Dfxuu5+Qzw+o44Z69ex4Q27fAZzHx84ed0iO5w3JreoJfnVj3H7F2ziQ4LTvd4cKoKq0wmIsJh/lVQxEi3p9HLCRCY/xiH1r6PL+wnNRjEVFOf3n5IrQ/EtZnF3PD+StrFWvhkXH9e3/wis3fN5vwO5/PY4MfI2ZdTofnjlhwmfraOcSk/sdG8gK0GjYsdLh4pLCrtByoTkQwj/lGlDm3PtnPvzA3szHVUfHZ2tyReuaYvUSZ9neU8lrYyctEbFAYy2XTroloPrb4gr7rmHe8/zO+67/ky30U3Zx2ZTGNSYVj92oo3EOIfczYzZ/1BAAJFB2kVreP5cWdzVtfU0peEsv4Wb3GtvyOSe7Fi8G08uvlX8kLrUJQwUYYoWllbEQgHyHJkEQgHiDZGc0WXK7i66/VM/mIPv+8u4MFBUZzTPfnYr20ddeho1PVcKS+nAJaaTXwdYSXDYMAgBAOEgSt7jqXzsAcarPnfdf/lvc3vcbvNxR1FhVR/HRPAtNho3o2J5rb2F3L3iOcarAnAmukU//wkS1UvOTod8aEQwz0eEkKHB31qaivHq7kiewXTfnucjZ7sKsfphGCE28OEEhvd/IGmeZYdrY8/hrbSHNS1IK8+zvETQB7wNVWzVRQ1oo1Ay3aOq1BD5xUGthkMbDAZydZphIGkUIhTvT56+vzoe18FV7zbIE3v9/fyUlwMs6Ii6e7z86+CQtL9gYpDgsC3EVZeioslrMDLCcMZcvFbDdJk0b8Rrjw+j4zgjdhobJqGIVxaZ/yqQmQozPgSGzfYHWgADSznwl8e5cH935AaCPJyXj6dAsEq33sVhddio/k4OooBHi//ySsgevDdMOr4QuPDIsyny5/lvZ0zKVJBEQKjEHjV0lGDzn4/dxfbGOH2VO1Mh9573JpAjXXID9hVFaMQRAhxROdd07W1uQNk5Dk4WOwhFBZYDBptYy10SY6osginXNP7/b18b7XwRVQE24xGANoEgrQOBlGALL2ObJ0OTQgudrqYkDyMtlc2YGfIGsrpVBR+iLCw3GwmW6dhDQv6en38zekiNRhscB0K7V/BpPnjWaEL82ZOHoNqcIwBhIBMnY4b2ySj16x8fPk82kQmH7OeO+Dmju+vZ51tFzeV2Lm9xIa1rF/dq9fxSGI82w0Gni4o5GKn+/CJDSync+XbvLPmZeaYNWxa6b0ufwCOK7HRs1LfUJvmtkN2rnlnObFWAx/d0ofn1j/K0oNLGd9nPJP6TjpyRDdrFfbPbyXStZ+gAu/FRPF2TDTJwRDP5hfSz1fDtS7TXLa7gAkz1mIxajz8f904tV0s87fl8MKPO+mXGssn4wZi0NUwWlcPZ7EyAhia2hGjNpBFY9+u9bi6slWwZjqZP/6di9q15rGCIq5yOI8ufJT7GQ4Lbv9kLfO35TJxZEcs8WuY88cs8rz7ATArZgb5vQy3FTDM7aVV6MiJ2kJV5bsIK7MjI9hn0GMOqXSN+RuPjbieLrFdKkY7A6EAK7JXMHf3XBbuX4hRM3JZp9Gs3tibTZmCt2/oz9ndk0tfludMgOI9tdpd7iFU1ISG9n01aBapKlMS4lhstRAXCtHf68OtKKw1GfErCmOCRu4751UMHeofQ16ZL+bfx7+yF3Kl3cE/C4trzWEbBp6Kj2N2VARPGNO44ppvjksPgKxVuOaM53WlmJlRkQQqtSW1rH+9t6iEhHAlJ7kBfUIwHOS5eTcxq3gjbQJBrrU7GOT1EhkOc1Cn41eLma8jrNg1jYucLu4rKiEpFKpT0xcMsbfAxb4CFy5fCE1VSIoy0i0lijhraeiO2L+Sfd9MYLMvj316HV5FQS8EbYIhevr8dPf7j7zeDfWHGomGOsd7a/hYCCGOunRZUZTzgf8CGvCeEKLOV7GW7BxXZKvI3Qibv6j43KEofBYdyZyICA7pS6NUzOEwCuAuc7DiQiEuczi5NBDN/FaPsjbcGU8gjEmnkhRlpGNCBJ2TIujeKpLECGPFA6ki6H3xw+za8TUPJcWz22BgrM3OPUUlR8TElDe9bE1jYkoie/V67nLF0PvCb+nbLga9dvRposqabP6CAPB0QhxzIiMY7PEwodhO37IH4GajgXdjolliMXO6x8vLeQVEh8MQmwaXv13vN9JyzbULxzPZs5OePj9v5uRVjERVr6EKMM9q4fHEeNoEgryRm0+71DPgxq9r1QiFBVlFbnbkOMjIdbBxdxbF/jxC5tfYo7czxO3heruDAV4fZiEoUlV+tZj5IDqKfQY9w90ensovrNqRHWc5y6+tH1hgtfCrxczaslH/cizhMJ38Afp7fQzyehng9WIUpZq+i99gvqM9c9cf5NeMfILhI9uwQVPp0zaaPokqAzvGc9qeZ/jqwEI+jY6kSNPo7Pcz2uHkbJeHlGoP4b16HV9ERvBlZAQKcEsJFFsfpDj2VNJTIhnQIY4+baJR1ZqnwquXszJLzCYeT4inUKfRJhCkfSCATVPZbjCgAFc6nNxXVIIlpuNxX9uXDs7no+ioKk5NTXWonC0GA7e0SiIxoOeM9h9x85B0kiJNR9XMyy9gTWYO722+nT1KCc/kFXKB+7DzW67pURTuTk5ktcnIG7n5DPVUGp07zjqU+/Nd3OtYT46mca7LzXCPF4MQbDIamBdhpUTT+JvDyf1FJcRXrrOVNHcbe3DNO8vRaypv3NiJZ9Y9yK7iXTw26DGu6HrFEZrx614tnXWoxgajgUcS4zmk03GjzcGtNjsx1TTt1vbcYruJwqTWTDjbSrRFR6eYTnSO6cx3G7O5Z+YGJo3szAPnpVfVrKEOHY0sncYF7dpwX4nglss/qPXaVneOq2uGgRGpbRjh9vBUQT3Hguq4nx8t28eUb7cy+f/SWO97heXZy+lm7Ea/xNOwZSzBGNzKcouJ7LJ+oE0gSE+fj9hwGK+ikGEwsMOgRygKfb0+rrQ7ONftxhiThlJHHdpr28u7m97lf3v/h6ZoaAV9Kcrty+J++2m77Z0qxwogV9PYbDSw2Whki9HANqMBn6KQHAzR3+vlMqeLftFd4IKXUdvXL2ylrjq0T6fjjpRE8jWNu4ptjLE7KJ9DKFFVpsVEMys6kj5eH70cFxHqNo4z0xMZ2DG+5pepaprr54/n776dDPN4+W9ufsVzs7Y+IQhMTE5kldnE2w6N0y9999j7oXWvsmfVNO5JSmS/XselThdX2Z2kBQIc0On4OtLK51GRxITCvJKXT1+f//CPHEOfIIRg+dZMVuzLZmP2/awzOLnRZufu4pLS50VZOcvLZlcVPoiOYkZUFBqCCSU2brA5MJRpiransSPHwW+78vltVwGr9xXhDYSP0A157HRIDHNe0vssUw+RpS+9Y1rZ4FJAUSpeBuJCIS5wurja7qRDsNJg1zH2fU1Bg5zjBohqQAYwCjgArAauFUJsq+2cE+Uc+4NhcmxeHL7SkZZIo55WMaYqzuTA//SE4r2sVEorgR+YFRXJOzFRlGgagz0eLnYHOc0fIsVrRwDFqsoak5HvIqwssZgxCsENNgdu39Usi7gUrz9Ets2D3Xu4wsRbDXRrFUm3lChS85fQd+fL7Eh08WJcDFYheCa/kCHlD1cBLmFAI4RJPezgKJQ67feWNfC7813kuEYT6Hsjl/RtQ7/UmFpj/DJ+/QrmP0ZX/SFsqsL9SaW/cXuxjTtLSveCEaqOsKJHC3kQwDcRVv6VEEfbQJA3c/NoEyyLBbzoP0edPhFC8NOcj/hjzfO8l67R0+fnrZw8IoSo0nn50QMCgzh8rdabjNydnIAKTMvJpzsWck+fTG6Xa8l3+MgsdLErz0lGbqlDXN7IFUXQWZlFbuRKLCl6Hi4o5hKXC4WqHaai6gioBmZaNP4bG41VCJ7OL2S4p9rUYy0jKUIIsm1eit2lHV/BgX0kLn6AVM9qPoqOYlZUBMWaRmIwyGleHx0DAaKESkBVyVYEOw16NpqMBBQFczjMII+XM9wehrk9/Md7M4sjLuSSvq0Z1CmedrEWDJqK3Rsgq8jN+qwSVu4twr59Jt2i/se6Libcqsowt4dbbHb6e30VHWZI0RFWdOjDh+uVokCOpvFkQhy/W8wM8HjpVzCAqc7rAUiOMnL1gHaMGdielOiqjmTlOlSOR1GYWjbr0cXv5/EiB6cEBErQjQDyNI33YqKYFRlB+0CQl/MK6BwIwkWvoAy4uc46VFlza6yNxxPjGWNz8EhRcdVRL81YWriQ/4jzfzWbuDs5kYHuAPGHLsDXeyyjB7Slf/tYjLrDI/Aef4g1mUX8sCWHbb9+TjDmC7I6mHgmv5ALXe7yy1dxHf1CA1VHUPEztlUyh3Q6PsnOIa3ajEh9R+PmLFrNgfl38ElXDzGhMC/kF9DX569STqfOxLtRFj6OsmIJh7m/qITLnK4qozcC+Ei5lGna9TxzdRwvbJiMw+9g6oipDGsz7Mhru+hZurK7yuei/H8UcGt6XoiP42uLAaMQDPZ4SfcHMAhBtk5jV5ljVz4jU07H6I48cvojfPmbme82HeKXv4+gXZyldNp27kS6ulfVfjF0ltL4VL+jysffWy08nJTAFwez6e4P1toPHRFWUUO9vTM5kUM6jbkHcw6fqOogXO3+Vafa/Sx2+TnjhUWckhqBtd3HrMhZwWODHqP3IYGy8HG66A5r/qHXs9xsYq3JyB96PcWahkkI2gcCDPB6OdflrphRO+z4HL2/zbJn8f6W95n1+yxCIkR6nEp3v5/4UIiAopCl05FhMFCgOzwL0c3vp6fPT2Q4zH6djuVmMw5NZYjbw6MFRSwSF7Omy72M7t+W4Z0Tan1hzvj1K1j8LF1F1Tp0UKdxY6tkgorCa3lF9AwoqCF3RSMqf0zNt5j5Z2I8UeEwIw724D33LUQadYzolsSoHsmMSE+sEpYTDIX5dslaChbeybud3XTyB/ggJw9L2XOlwkrNCIoKQU8Vu2yqwnWtUrBrKp8fzCGy7+1EXfxMrde2ejnD+kzGpyShAC/lFVQJ6xKaEVVRyVCC3JOcQK5Ox1P5hVzkclf9sTr6hBK3n5mrs/hs5X7++GMz3Vu9SV4HjclFxVxnLxsQqPQwCygaIXSYyib/D+g0XoqL5RerhdRAgIcKixnm8fG8OoF33GcA0CUpgqGdEzg1NYa0hAiizDr8wTDLDqzhy8VT+SOwE10rEwPdXka53Zzm9dI2EMRA6Qh8rqaxzmzhF7OBX6wWgorC2S43dxTbSA9UmtUaei/inCdOyCLIho4cW4D7gVQhxHhFUbpQup30vKOcNxh4QghxXtm/HwEQQjxb2zknwjm+dfpqftmZR/XLoKkKHeItDGgfx1j3h9yU8SoAK3QR/GC18FpsDAf1OgZ5PNznDNNjxOMw4CZybF5W/Po/eq5/is7hskF3BfbrdLwWG81PEVZiQiHGmTtyzVVzMKgGilx+MnKd7MixsyPbwY4cOxfmvkUP23e8GRvNprRohro9PJ1fSHwoDArYtDg2d7kTy+Bb6dk6GtPGGbDoGXDlVpTBp8CDiQksslqYVFRCUmEv7gtOom2smfN7pnB+rxROTY1FK+/QFkwh49upAJiSjUxMTuSgXscTZdPBbkx8EhzFs8FrAeinZPCU9iE9tUzWmI3ck5SIWYR5JyePtEDpNH1NcWqhsGBtZjE/bc2hw/rnsZR8x5SEOHrGq7yfXeoYV2CKhXOeqOj0/Svfh0XPovfmA6WjnHcmJ1GoqTyfV8hIt4fFod7cHHwEgIQIA+kpkaQnR9EtJZKkWD9frr6dBfsz6eP18armoHWw9MWiYqQvbGB66Fz+l3I7F/dpzeWJh7D99hAPqUXsqjRyf7grBpJ74zjnedaLrmzMKmHjgRI2ZNkocJZ2RteoP3NJyQyWR6l80zWeIk1jpMvNtXYHp3t9aNXKSdYqQt/djy9/C2tNRpZYzPxqMVeMKvX0+Rge2Ynh5zxP+6j2RBmiUBQFV8DFPts+1uSuYdGm6SzbfxBNCP4WGeIWm510f6CirjuFiRmhc3ghNAaAcZYljFfmkhjKqVwy5kZYeTY+FoMQPO63Yug8hY+yklickY+qKJzfM4WbhnZgQPtYlIVPVNShrvGlD9mtBgMPJ8aTqdcx1hXgrtMeRO1/C6v2FrF55QLO2PUs3UQmACtNRh5JSsCnKLySm89pPh8rLGezZdCLnN4xnp6to46cAVkwhZ3fTmWzwcC/e7fmNK+XN3Ly0Sh7CFa/tmumH9FWAL6IjOBfCXFcZXdwRn5bxvofwaRXaRtrIcKow+4JkFnkJhQWPGT4jMXab2w0GnlJc3Jx9QebzoT9lFuZGh7Dpyv3M9y8h8mRHzEhykVUKMxnh3KOjNGtIwZwV66DRZ+9QNy+d/h3ciSpiTrezckjKRSqcFKLRARfxtxC5/MncVa3JPYtf5mntk9nrR5O9Xp5qLCEXv7SFwMBuBWFD9r15gO9kzhTHK+f/Trd4rpVFV4znYyPJlW5n+UIwKcY+cA/CnHOk9w5ohN/7JjLrJXP81vQxiGdhlBKQ6+6BPx09/np6ffTJaID6pkPskUHH275kH32fUzs8yCvzE7koj6tmZq2nozP/wF++xGa1HQ/s1bB9/dDTul2z8/ExTI30sqyzAOHZ9dqmMat4hxX6vsqa74ZE8WbMdEsy3MSkTr08KKi8kVktv1H2ldOpfv57A/beWfJHsZesIPZe6bz5JAnuXzfxho1RZn3JnQm1IG3VzhIYv9Kir+8i1j7zgrv7gh34mhxwQumsPK7V1hqMbOnfSS79HocqoqKoE0wSOeQSs+08+jdbxzpcekY1n+Gf+HTpf2tKH3JnR1VGmIH8GR+IRZ3V673TqZP22gev6gHAzrEVdWspQ7laypjWyVToml8mDyKwvR/8uCXm0i2b+SF6Dl08myqUr6dBj13JifiVFWedJn5I+p+PshMpNDlR1WgTayZOKsRXyDEgIJv6Of4hGeSI0hM0jPjUG7VWb+j1CGAPXod17VOoW0gyEeHcsmKPA31+jl0ST5yRaMQgqLimunQAAAgAElEQVQl75L/1f1k6A0816sVVhHmvew8UsteZNzCwEehc3nbMJZhnRO4MvkQvfa9yoPBvawxm7i/qJibbI4617jszHEwfdlevl5/EG8gzL0pv/K9czYHdCova07+z1V15iqks7K1zVV8YB7LxqwS4oo2VDyzAZaZTTwXH8s+g55hbg/3FpUQE9EH3blPEt99eMVv2Xw2FmQuYOaOmews3onuoIez3W4mmn10LH9RK+vO3KoFZ6shxIx6AGPHwbBmOgW/T2UWNj6JjsSpqpzndHFHia1ikMBuTiV6zPvNPorcUOd4FrAWuFEI0UtRFDOwXAhRZ44cRVFGA+cLIcaV/fsGYKAQYlK148YD4wFSU1P7Z2Zm1rNYjcPqz56k4/7Z6HUaWV1v4kDa1dg9AfYXuUnc+Tn/V/QBicLGwKATh6rSs10aO40Guvn83FNcQrsON7Kq091k5DpY9kchWw/ZAeiXGsNtHQsZteNxdLbDkSlbDXpejY1hmcVMcjDEDXF9ufD8V0kwJ5QekLWKgrnj+SpUyFseIyHgIaOXMXZHxciPUtcoU7X4vADweGI88yKsjC2xM9EJsyJv5pncgfhDYawGjbuil3KN+xOiQ0VkFIbYajDwSs/SXaFeyS2gn9fHB8ol/Nz6Tnq1iSI1zkKkSY9A4PAGYf8qRmVMwa7mc3tKEgEF3szJr/IgzrSewneJE/jV05Fth+xcEprPg7qZ5Jj8XKfEEBcK81W4mLjKHVg9y1mgqtyVksg2g4GbbHZuL7FjUI0EBozH9H//KrVBCOYufICpB37Ci+DqjAIudrroVv3h2/sq8s59je82ZvPtxkNszCoB4JR2MVwYt5cDJc/ztQV6eX28kF9Au2CowqkWAt4KXsSL4TF0SoygT9toRouF9N3zJitVF08GLBzS6xgaLbi/qIReZSN+bwUvYmb0ONrHW4g06fAFwuQ6vOzKddIztIOXdG/RUctBALv1epZYzCy2mNlkNCDK3rYNamn8lz98eFS0k9/PgD02znJ7GFKWcUuU/e1Iu5l1Xe/DoKm0iTXTPt5Cmxhz6dt71ir4egIUHY4J3KfT8VBSAtuNBq62O3gg/Xry+v+DT1ZmMnPVfi7w/8QDhjnEU0xGYQgF6BSv8X50FG/GRhMfCvFk4pn4059i/tZc5m/NodDlx6zXOKtbEpcmHGDYln9icmRySKcxMTmRTL2+YkTFHjazPNyTj9RL0FIH0i7OQmffVi764ymSggdZXKzwQHI8bRL1zMjOIbo81OQY2grAy7ExfBgTxT1FJdzk9LMicTSfWG/B5Q8SZdLzt9B8Tj/0Ln+PUVlaAvcU27jTXC3Otoa43smzN7H5oI0ru2xkgfY5Z7jd/CevoObYx0oOjjcQYsnnL9Jv96sE9F4uV0tv5JxwcZWYVM9pdzEj8mamL93HIZuXTolWbhuexgW9U/h50zSmbvsQm1paJzr7A7hUlXUmI25V5Vy3l392vILY8ytFvVVyFjIKS3WOcFR7X4Xvkrd46KtNfLPhEIPT4pk4sjOpcRb2bVxE6yX3007koFeouZxD78Uz8mEmL5nMoqxF3BAezg175pGi2dhVm2Zd97Os3l5tcmENCz7Iyav6vc4ElZzNEdNHgM/BYl8YivbUWM7fzCbuTEnivXPfY2CrGkIIamgr1fF1H80pW0YzpIeHtYEnuSi6O//OWAeu3Bo1w8C3oaG8YPk7d53dhcFp8eQ7fbz32x5+2prLzal5PBp4FV0dccJHOMmV7KysWXn039n1CiLHfACA2x/klx15zFieycq9RYw1Lubvuq+IDBWhUDri+2BiAptNRm4tsXGHzcdMzmeK+yqu6NeWf1zQjfgdn1e8QFQvZ4mqcnOrJA4aTEwbNZ3/rTEwfdk+OiZYeenKPvRvH1ejw5qnaUxKTmSnQc/kwmKuie3Dzl5/5wdbKnsLXJyaN5dLbR9RojkYQwwAX4oS2lWezq9HHSq/n0vMJiYlJ3Key83zeYV4hIG5hotYlnYXiRFGQmGBOXctV2RPpavI5GuHxpOJcbRK1PN+Tm7FwIur6xXsHv4yu/OcrNhTyG+78sm1+9CpCmPaHsSmvsovZrjeZufBopIjZnl2d7mVx11XsXxPIUadyrPt19LH9hGTYnQcyg/wj4Jirok8+mxUscvP1kN2wlkr6bP6YaI9+wkAn0dF8mZsNC5VJc0foLfPR0xECo7Wp/JH0M6Wgi2ERIiuIbimuJCu+xyYEVVf6oClKddzT/6lFLr86FSFrsmRpESbMOpUWjk2c0X+s/wc5eHT6Eg8isK5Lje3ldjpGgigqDq4+YdmdZAb6hyvEUIMUBRlvRDi1LLPNgohTjnKeVcC51Vzjk8XQtxV2znNPnK8Zjr/Wv4ExZrGBU4XwzxeTKqhNMVKOAxBNwWqymKrmfHOIryKypmtU5lYZCPdEcnk4B2sE6WrNfWaQr/UWIZ3SeDCPq3pmGA9rDP7tiPi5labjEyLjWadyYQiBB1CYeJDYYoVwR59aVzZKXtKuNXmYGRM2T06lhidjy+DPb8ApZ3tM/GxzIqK5Aq7k8cKi1BUAwFFjxLyoedwo/qv18Q7MVGkxWu8nptPvLENzv+bRnLP4Uef9ph9G1nbZzM+JYkiTePV3HwGen0VHbBQwIcBDYGeABuNBiYlJyKyfbyQV8DQ8pSpx7KKtqycbkXh+fhY5kRG0DoQ5EqHk74+H6qisdVkYY5Fz26Dnn5eL1MKigjmlDo0hxt3zdOSmYUu5m3KZtGOPDYftOELhrk15kV+SMpHKDCx2MbVZTFyFU6yZkTV9BAKsFEneC02hpVmEwlZLm4psXN9RAAFCMR0ZF2/Z1kV6Mz2HDuHSrw4vAFMeo04q4FuZfG9w7skYPnujiPqUJGqstZk5JBOR6FOj6LpiQyFSfU66ePzkxIKHfnwPZY6VK3eBoD/xsXwUXQUnf1+ni5y0RM9YZ8TpcwpV4CMwhB7dHo+6p7IFqORke4wxsCdfJuXhi8YxmLQOLt7Mhf0SmFEehJmg3aEpl1VuK8spOe2EhsTi22oZdfYh4EwYBF+KAv/uFqJwa0qfB0uexA2oK08nBjPDxHWwzMEqg50Zgj5yFLC3J+UQIZBz53bchnp8dTr2gZDYd5esof/LMwgImk5wZi53FVUwnibvVazgqqRQDiMSQSwawo3tUpmb0GI5/IKOT86XKNmIBTmf5uzefvXPWzLtqOpCu1izQRxM4TncERks1+vwygEp3h9XOZ00ac83rG8nOEABA+HDh3pwFVtK+Gw4PPV+3nhx53YPIenS9MSrHwa/wGtMutY1KQZ8Wt6xiZEslev8cWBHFJDweOut56gh0GfnM4tNhv3FNtq1UTTMyKQDyLMYqw1lzM2jfda38d/3S9y5yl3cUff8bUL19DHlyMAX1hjXNskDupUvss6VDE7Vls5NyldeWTO5orBFgCLQWPiyM7cfman0tm+OjQr0JkAFYKHRxSrO8cOXTzPeS/js+BZtIkxY9SpHCj24A+FaRtr5sbB7bnm9NTSsIVKmgFKnytfRUVyptvDs3kFmIWKO2zAoPoxlmV+VappOhSFca2S2W0yc98pL/H+Ah17C1zcNKQDk8/vVrU/gNIsCJXilN2KwuTEeBZbLYyxObi3uASzZgRURNDNLxYzjyXEE8jx8kx+EaNqaSt1Uqmc70dH8p+4WCZVaq8BNJyY0REkktK28pvZxMRwJAmhMJ+LkrK1HDU/V8JhwaaDNn7cksNPW3PYW+BgVMozrIh1cq7TxTMFhYdjhsv+61N0KKoRHQE26+Ge5EQEMHlLDl0C/gb38YWqyo8RpWtg/iibUTALQYdAiP7+ICMd9orsTXW1z0AozO+7Cli1r4jt2XbyHT78wTBmg0braDN3Fr9A25L5fBIdyWdRkbhUlXuKShhns8PZj8Pwvx/d7kaioc7xMuBsYKkQop+iKJ2Az4UQdV75kyKsYsZlvFq0ltmRERRpGpoQpAUCxJWlWMnWaewvCzTv7Pczxu7gUocL0eNK9p7xCsVuP5qqkBhhpE2sue4Fb9UaeDl79DrmWy1sMxiwqypR4TDd/f4qcWXA8a3urKQpgNfK0tT093h5tLCYLpXifrI1jRfjY1lgtTDM7eG5/AKie155XJr5K15lfEoSmXo9DxYWc43DWWWqSAA/WC1MSYgjKRTirZw82pW9YR9XOSuNAC43GXmr7KWjMuk+P2Ntdi50uY8cxapniqJAKMz+IjdFLj+hrW/zftb7rDCbaRUMconDRX+vl/hQGLuqstFk4GeLhU0mI9GhEHeU2LjK7jwcipF2Vp0LCGukljpUb4732s67p8pHS80mHi1bVDfS5eYSp4tufj+qgAyDnu8jrPxktRATDjOqMJm5vkfpEG+hV5toRnZLYmDHuCMzaVSmrJwB4N8JccyOjGCEy82UgqKqU6OULqa7PzkBp6Lydk4evf3+BtehIPBc2ctkR3+AqxxOEkIh1hmNzI6MwCgEz+UXcEbl2PN6am47ZOe+L9aTqb6HPmoDr+fmVf2dMiqP6PkVuD05iQ0mI2/l5DGwPH6xjjokhGBNZjG/7sxnf1GpY9Q1OYJrtF9IWP5svTM/HEEdL64ef4jlewoocPrpEG+lX2oMOk2tsQ5VJ1vTuKxtK07x+ngrN7/q1PIxtJVlB1cxYeGtTPL3ZMLBH+o8dgSlKafKneMqlN3PHzZn88DyMZzWphvTL3yzbvFayimABRYzf09O5In8Qq5w1pLqqlodEkKw8YCNjFwH0WY9g9LiiTbrq57TSH1Cjs3LNxsOsiPHUeEUD++cyOBO8YfD7iqXs6ytCGBWZATPxcfSPhDktdz8ilHamoZSXIrC7e06sEUHvfX38tumRFpHm3lxdB+GdE6o3c5q1zYEvBwXw8fRUcSFQpzjchMZDrPaZGKTyUi6z89/8/LL1r9wfP1tWTmFt5iHE+P5X4SVO4tLuL3EfsTzbFZkBM/Hx9LFH+DNnLzSRbD1HOQRQrArz8mPW3LYuvMxllm3lmViyq+YAavspc2JsPJMfBzJoSBv5ORXXeR2sjxXykbobSX7+CQqkuEeD30C4ZNu5Phc4FGgBzAfGArcLISoPYFk6Xk6ShfknQ0cpHRB3hghxNbazjkRI8fMu4cgsMpkYo3ZyE6DAZuqIoCUYJBu/gBDPR5+djpRUJh0+ZvHn6OvHlNwlSl0h8EYQ/wl/2o0zbkRVl6Ii8WhqfT2+kgNBsnRNNabjOiF4LqDxVxvd5J41asN0rR9fy+TRR5LLWZ6+nxcbXfS2R/gkE5jdlQEy83m0pjf3HxwBhteTqgyApinaWQY9ChAmj9wRFqkQncYUIhvQDnF/pX8/uPdfBwuZIXZfMT36T4/lzmdXOZwYRGiTBPiR91//KmQTlQdqja9aVcVPo2K5NOoyIo0YuUY7X4ucbi4+6xHiR50x/Frfj0BUbSHz6IimBoXiyksGO1wcrrXiwB+tlr4JsJKYijEE3tzSVesjVqHlphNvBIXw25DaciKXgj+z+nirmIbKaHQcdchXzDEqz9v46PMB9B0BUzLCTLMv6/i+8oP3gBwf1ICi60Wns8rYGC+s9HbytFotHpbrQ5V59OoCJ6Lj+Ol3Hz6F5QukD1WzX/9No0v9rzNlD5fMjrBXmdbqe4c19RW9he6OffT24mJ28fy6349+uxZWS56MpcBh52aq1un4FEUvj6YXSXLUGP0Q/W5tpVplD4BqtShlSYjf09KQFC6AG1w5QVoAgo9YXI0jQe6dOOg6sFz8FpUzymMHdyeu87uUmUxXa3UUM51RiMfRkey2mzCqyh08gcY7XAy2uHE3hj1tqycgT2/MCUhnu8irfT3eBlvs9MhEOAPvZ7p0VGsMpsY5vbw8L7SNTMN0fx+9av8c+u7RIVD3FJiZ7jHgyUs2Go0MCMqkrVmEwM9Xl7MKyA2HG629lmZRqtDC6bApi8grkPTbkxSCw3OVqEoSjwwiNJ+e4UQoqCewhcA/6E0ldsHQoh/13X8CclWUc8cmgNVBaLasPLeWn37RtfMiD8Pznu6aqLt46VSgu5iVWV2ZAS/WUzkaDriwiGGuL2Mdjhx6HrBkLvpeuboBkuK1R8yd8XzvGOGA/rDnV9iMMhYm4Pr7A50OhMZba6EoXc3Tjnre22NpzRaOVkzHdvPT7BNuLFpGpHhMJ39AZIrO+SGSDLMp0K/sY2m2ex1qAbNALDRZCRTp0Mo0C4QxCrSMQy5p/HKuegZ9vgLmRYbw88WM+HyOOuw4HKnk4nFNvJiRzVJOctTW9lVlbbBIJZK/WVD69C2vD3c+NN1iJCZC7J7crfvWxIoqfjeqyg8khjPz1YL/ywo4mpP8MS0FTrDyEeavN6GgCvbpOBHYcI2C+mX/uOYNS/7ahwZxbtZeNUPtIo216lZ3Tmuqa0IIejz8uOQMJf5V8ynVUSr+hlS6eVurdHIza2Ta8yX3Nj9UE0LTavTVH1Clk7j7uTSFKJXOpxcZXfSKRCgUFV5wd+WedYQuhQLqcEJXNB5JJf3a3P4Hh2rZg3lrJKFgsa/tmLhE8zR+XklLqbKgEB8MMSEEhtXO5zsbqS2srVgKy8ueoC17gNVPk8KBhlfYudKh7NiBvSkr0MnkIaOHP8shDj7aJ81Bic0z3H1bTOrbac7cF5pqHSVHfIaSvXte6tthZrhKV1Q0KgVsPpWndW2CG0KzfDqD9m7+nWy/CXEBvz0CITR603Q9zoY9WTNu+00lBNQziPqEAo0Zzmbqw5V37632hahTVZvV7yBzZ1PhhpGQaFbSBCRdHK3lU35m7hj4R3oVB0PnfYQ5xfmoq58k7X+Ip6NMrBLpzHZGeC69Gv+3G2lrL/9OTGVe4OZ2DeN5NMx4xjWv0+9fzIswgz4eCghVw823PHekaO81ba1HREsABQWtxtZZx362zsz2Wv8N88Pf54L0i6otz3eQIhn/zWZglY/sMEEC3NsmEUzXNvq2/dW25K5KduKy1vE1GgrX1sMBBUFHaV5g305Pk5NOpWXrniJdlHtGk+z+nbMlZ7ZTVVO729TWR12kKfXkxIS9PcFMDXRtd1n28fG1a/jzfiRNLedvv4Q+uZonyeiDp0Ajss5VhTFBFiARcAIDr+URQE/CCG6N7ahLXkTkBq3j25imuRBKDWlptRsUZp7Svbwj9//wdbCrUToIzBoBoq8RSRZkpgyeApntD2j0TWPhebUFEJw6dzRbN26nzHt/8VTN5xf73N3FO3gyu+uJE2M45ub6o5xhqNvH13Ok99t5sv8WxidfiFPDH2i3vZ8vymbSV8uJLLzS9zc6ybu63/fEcf8We9nriuXpYeWst++n2hjNG08begQ3eFPV06p+ed1jqtvslaZCcC9QGtKU7lVbLICvN6oFkokEslflLSYND654BMWZS1iZfZKguEgvRJ6cUHHC7DoLSfavGZFURTuG3AX4zPG879dP/MU9XeOf8taDsDAVo0bt9i7TSwzMzuw7FAdG5LUwDcbDhKdvBoUhWu7XduoNrV0kq3JXN7l8op/lztTEsnJQn3CKu4SQrzWHMa05JFjiUQikTQ9QgjOnvk3chwuvrx4Nj1bxR79JOC678azIXsX086YxchuSUc9vvrIcW3szHFw0cf/xpT8PQtHLyTZmnzU37Z5Apz27++J6PosZ7UfzktnvlSfIkgkkmakrpHjOnKPlSKEeE1RlCGKooxRFOXG8r/GN7NlM+33aUz7fVqzahYWFlbs1S41pabUlJp/BU1FUbiuw9WIQDZvrKwjR3IlAuEAO4o3EnSn0adt9HHp1lbOTolWNF9nAFbl1G/0+KctOYjINQSEm+u7X3/Mmk2J1JSaJ6Nmc3NU51hRlBnAS8Aw4LSyvxo97T8zMzbPYMbmGc2q+Vep9FJTakpNqVmZIbH90Rwx/F4wi3C13NY1sTZ3LX7hJlbpQ3yE8bg0ayunTlPpHp+OKiyszlldr9+auzELS8JyesX34pTE2vfL+qvcT6kpNU826oo5LmcA0EPUJ+ebRCKRSCQNRFM1hqVcxFJlJh9vnM9Np9Yde7xo/yII6xjUalCT2DOwYyI7/+jAiuwVCCHqzHecZ/eyKmcF5nZ5XN/j/qPnRpZIJC2Oo44cA1uAlKY2RCKRSCSScm499WJEIJrpW9+v8zghBAsyFxF0dWFoWusmsWVQWhwBZzeyXdlkFNe9uOy7TdnoY38nzpjAue3PbRJ7JBJJ01If5zgB2KYoyk+Konxb/tfUhkkkEonkr0u02URn40UUhnawMrv2hdo7inaQ78km6OzOoLT4JrFlQIc4wq4egMLP+3+u89hZG1ehi9jFdT2uRa/VY+c3iUTS4qhPtooza/pcCPFrYxsjs1VIJBKJpJx5m/YxefU19EroyaxLah5Bfn7V83y6fSbW7CdZNvnieocx1DdbRTmXTPudbMtU2iWozP7b7BqP2ZXr4OKZE7HEbmfR1QuJNh7f4kCJRNL0NDRbxa+V/yjd8OaqxjaypSOzVUhNqSk1pWbzao7q0Q7VfgbbSlaxvXD7Ecd5g17m7ZmHcPZgaKf2DYrvPVo5h3ZOwFbYjYziDPbY9tR4zIw1G9BFb+RvnS6rl2P8V7ufUlNqnizUJ6wCRVH6KorygqIo+4CngSN7qT85MluF1JSaUlNqNq+mUacxqu3liJCJF1dPpfpM55xdcyjxleAuHMQ53Y+e27g+mrVxbs8U/CWnoKIxJ2POEd+Hw4J5+z5HQeH2vrc0imZTIDWl5smo2dzU6hwritJVUZTHFUXZDkwDsigNwxgphGjeIVSJRCKR/CW5qn8XfHnnszp3JV9mfFnxuc1n493N7xKvS0fzdWJ4l8QmtaNPm2iSrYnEcCpf7/4ah99R5fvPN6zGb13KaQnnkmKVa9glkpOZukaOdwBnAxcLIYaV7ZIXah6zJBKJRCKBgR3jSDOdjTHQjedWPceCzAW4Ai4e/u1hSrwl2A9cxIj0JKzG+mQmPX5UVeHcnsnkZg3B7rczfev0iu/CIsy0jS+hCAPPjnywSe2QSCRNT13O8RVADrBIUZR3FUU5G5AJGyUSiUTSbCiKwq1DO1Gw5ypamTty/+L7GfzZYJYeXMrl7e+iqDiR0f3bNostl/RtjdfZmh5RZ/LB5g9YnbMaIQSPL3kBp7qd4fFjSbY27Qi2RCJpeuqTrcIKXApcC5wFfAR8LYSY39jGyGwVEolEIqmONxBi5EuLSYjUGDuqgEOuQ4xoN4IpXzrItXlZ/OBIDLp6LaGp4FizVUBpTuWLp/2OM+DA2uF1DjkPkWxN5qDzINgH8evN04g7zh36JBJJ89LQbBUuIcSnQoiLgLbABuDhRraxxSOzVUhNqSk1peaJ0TTpNR46P53NB1zkZffl7n53k7E/jvX7S7hzZOdjdozro1kTiqJw11ld2JcnOCf6GcZ0H4NVaYfn0Ggm9Zl8zI5xS7i2UlNqngyazc0xBWkJIYqAt8v+/lKUZ6qYNGxSs2mWV774+KZJbC81pabUlJoni+alfdvw8/Y8XvhxJyv3FLF8TyGD0uK45rR2TaZZE+f2SOa8nsn8d8EBzuh6Jut39WBQWjy3DktrMs3GRGpKzZNRs7lp+Ou2RCKRSCRNjKIoTL3qFCackcYf+U7O7ZHMm9f1R6c172NMURT+e82pXHt6Kn/kObnm9FTeuXFAs9shkUiajqZd3iuRSCQSSSNh1Gk8ckF3Hrmg+wm1w6TXePby3ifUBolE0nTIV12JRCKRSCQSiaSMo2araE4URckHMk+AdAJQcAJ0JXUj70vLQ96Tlom8Ly0PeU9aJvK+tDxO1D1pL4SoMfdii3KOTxSKoqypLZ2H5MQh70vLQ96Tlom8Ly0PeU9aJvK+tDxa4j2RYRUSiUQikUgkEkkZ0jmWSCQSiUQikUjKkM5xKe+caAMkNSLvS8tD3pOWibwvLQ95T1om8r60PFrcPZExxxKJRCKRSCQSSRly5FgikUgkEolEIilDOscSiUQikUgkEkkZ0jmWSCQSiUQikUjKkM6xRCKRSCQSiURSxnE5x4qifKAoSp6iKFtq+V5R/r+9u4+5s77rOP7+UNg0OFeFbiEFB24NDxrpsEEMiekedIU/VpehKVG3GzHVBNxm1ARMdLj4x/CJOIMYtrGyRQcEtli1GUw2JMYMWrTymLoOp3TUwYDhcAtL4esf59dyvLlPOb3u+zz0nPcrOTnnejzf9svv4tvfua7fL/lIkr1J7k9yzvLClCRJkkava8/xNmDTYbZfAKxrr63AdR2/R5IkSRqbTsVxVd0NPH2YXTYDn6yeLwGrk5zU5bskSZKkcTl2ROddCzzWt7yvrdu/eMckW+n1LnP88cf/+BlnnDGikJZnz1N7ADj9hNM7HNw7ltM7HKvxM1/SzFvWNV3SUe++++77RlWtWWrbqIrjLLFuydlGqup62uwoGzZsqF27do0opOXZuG0jAHct3NXh4N6x3NXhWI2f+ZJm3rKu6ZKOekn+c9C2UY1WsQ84pW/5ZODxEX2XJEmStCJGVRxvB97TRq04D3i2ql52S4UkSZI0TTrdVpHk08BG4MQk+4APAscBVNVfAjuAC4G9wLeBS1Yi2Ela1k9v/jwvSZJ0VOhUHFfVxa+wvYDLOkU0pbbt3gbAwvqFDgf3jmWhw7GSJEkaG2fIG9K23dsOFchHfvC2lwpkSZIkTS2LY0mSJKmxOJYkSZIai2NJkiSpGdUkIDPH0SokSZJmnz3HQ/KBPEmSpNlncTwki2NJkqTZZ3EsSZIkNRbHkiRJUmNxLEmSJDWOVjEkR6uQJEmaffYcD8kH8iRJkmafxfGQLI4lSZJmn8WxJEmS1FgcS5IkSY3FsSRJktQ4WsWQHK1CkiRp9tlzPCQfyJMkSZp9FsdDsjiWJEmafRbHkiRJUmNxLEmSJDUWx5IkSVLjaBVDcrQKSRAyR4sAAAlrSURBVJKk2dep5zjJpiR7kuxNcsUS2xeSPJlkd3v9yvJDnSwfyJMkSZp9R9xznGQVcC3w08A+YGeS7VX18KJdb66qy1cgxqlwsDBeWL/Q4eDesSx0OFaSJElj06Xn+Fxgb1U9WlXfBW4CNq9sWJIkSdL4dSmO1wKP9S3va+sWe3eS+5PcmuSUQSdLsjXJriS7nnzyyQ7hSJIkSSujS3GcJdbVouW/BU6tqh8D/gG4cdDJqur6qtpQVRvWrFnTIRxJkiRpZXQZrWIf0N8TfDLweP8OVfVU3+JHgas7fM9UcbQKSZKk2del53gnsC7JaUleBWwBtvfvkOSkvsV3Ao90D3E6OFqFJEnS7DvinuOqOpDkcuB2YBVwQ1U9lORDwK6q2g68L8k7gQPA08DCCsY8EY5WIUmSNPs6TQJSVTuAHYvW/V7f5yuBK5cXmiRJkjReTh8tSZIkNRbHkiRJUtPptop55GgVkiRJs8+e4yE5WoUkSdLsszgeksWxJEnS7LM4liRJkhqLY0mSJKmxOJYkSZIaR6sYkqNVSJIkzT57jofkA3mSJEmzz+J4SBbHkiRJs8/iWJIkSWosjiVJkqTG4liSJElqHK1iSI5WIUmSNPvsOR6SD+RJkiTNPovjIVkcS5IkzT6LY0mSJKmxOJYkSZIai2NJkiSpcbSKITlahSRJ0uyz53hIPpAnSZI0+yyOh2RxLEmSNPs6FcdJNiXZk2RvkiuW2P7qJDe37fckOXW5gUqSJEmjdsTFcZJVwLXABcBZwMVJzlq026XAM1X1JuAa4OrlBipJkiSNWpee43OBvVX1aFV9F7gJ2Lxon83Aje3zrcDbkqR7mJIkSdLodRmtYi3wWN/yPuAnBu1TVQeSPAucAHxj8cmSbAW2tsXnkuzpENNyncgSsS0llyyjxvffB0dq6LyMhPlaymRzokHMS0fLuqYfnjmZTuZl+kwqJ28YtKFLcbzUlaQ67NNbWXU9cH2HOFZMkl1VtWGSMejlzMv0MSfTybxMH3MynczL9JnGnHS5rWIfcErf8snA44P2SXIs8Frg6S4BSpIkSePSpTjeCaxLclqSVwFbgO2L9tkOvLd9vgj4QlUt2XMsSZIkTYsjvq2i3UN8OXA7sAq4oaoeSvIhYFdVbQc+DnwqyV56PcZbVjLoEZjobR0ayLxMH3MynczL9DEn08m8TJ+py0ns0JUkSZJ6nCFPkiRJaiyOJUmSpGbui+NXmgpb45Hkq0keSLI7ya627geTfD7Jl9v7D0w6zlmX5IYkTyR5sG/dknlIz0da27k/yTmTi3x2DcjJVUm+1trL7iQX9m27suVkT5J3TCbq2ZbklCRfTPJIkoeSvL+tt61M0GHyYnuZoCTfk+TeJP/W8vL7bf1pSe5p7eXmNsgDSV7dlve27aeOO+a5Lo6HnApb4/OWqlrfN97hFcCdVbUOuLMta7S2AZsWrRuUhwuAde21FbhuTDHOm228PCcA17T2sr6qdgC069cW4EfaMX/RrnNaWQeA36yqM4HzgMva371tZbIG5QVsL5P0PPDWqjobWA9sSnIecDW9vKwDngEubftfCjxTVW8Crmn7jdVcF8cMNxW2Jqd/GvIbgZ+dYCxzoaru5uVjkg/Kw2bgk9XzJWB1kpPGE+n8GJCTQTYDN1XV81X1H8Beetc5raCq2l9V/9I+fwt4hN7MsLaVCTpMXgaxvYxB++/+ubZ4XHsV8Fbg1rZ+cXs52I5uBd6WjHfK2nkvjpeaCvtwDUmjU8AdSe5rU4oDvL6q9kPvoge8bmLRzbdBebD9TNbl7Sf6G/puOTInY9Z+8n0zcA+2lamxKC9ge5moJKuS7AaeAD4PfAX4ZlUdaLv0/90fykvb/ixwwjjjnffieOhprjVy51fVOfR+frwsyU9NOiC9ItvP5FwHvJHeT5T7gT9p683JGCX5PuA24ANV9T+H23WJdeZlRJbIi+1lwqrqhapaT29W5XOBM5farb1PPC/zXhwPMxW2xqCqHm/vTwCfpdd4vn7wp8f2/sTkIpxrg/Jg+5mQqvp6+5/Ni8BHeemnYHMyJkmOo1eA/VVVfaattq1M2FJ5sb1Mj6r6JnAXvXvCVyc5OBld/9/9oby07a9l+FvLVsS8F8fDTIWtEUtyfJLXHPwM/AzwIP9/GvL3An8zmQjn3qA8bAfe057EPw949uBPyhqtRfervotee4FeTra0p71Po/cA2L3jjm/WtfsfPw48UlV/2rfJtjJBg/Jie5msJGuSrG6fvxd4O737wb8IXNR2W9xeDraji4Av1JhnrDvi6aNnyaCpsCcc1jx6PfDZdr/9scBfV9XnkuwEbklyKfBfwM9NMMa5kOTTwEbgxCT7gA8CH2bpPOwALqT3EMu3gUvGHvAcGJCTjUnW0/up8avArwJU1UNJbgEepvfk/mVV9cIk4p5x5wO/BDzQ7qME+B1sK5M2KC8X214m6iTgxjYSyDHALVX1d0keBm5K8gfAv9L7hw3t/VNJ9tLrMd4y7oCdPlqSJElq5v22CkmSJOkQi2NJkiSpsTiWJEmSGotjSZIkqbE4liRJkhqLY0mSJKmxOJakEUtyQpLd7fXfSb7Wt/zPI/i+hSRPJvlYx+P/qMX5WysdmyRNu7meBESSxqGqngLWAyS5Cniuqv54xF97c1Vd3uXAqvrtJP+70gFJ0tHAnmNJmqAkz7X3jUn+McktSf49yYeT/EKSe5M8kOSNbb81SW5LsrO9zh/iOxaSfCbJ55J8OckftvWrkmxL8mD7jt8Y7Z9WkqafPceSND3OBs6kN2Xqo8DHqurcJO8Hfh34APBnwDVV9U9Jfgi4vR3zStYDbwaeB/Yk+XPgdcDaqvpRgCSrV/oPJElHG4tjSZoeO6tqP0CSrwB3tPUPAG9pn98OnJXk4DHfn+Q1VfWtVzj3nVX1bDv3w8AbgIeAH26F8t/3fZ8kzS2LY0maHs/3fX6xb/lFXrpeHwP8ZFV9ZxnnfgE4tqqeSXI28A7gMuDngV8+4qglaYZ4z7EkHV3uAA49aJdkfdcTJTkROKaqbgN+Fzhn+eFJ0tHNnmNJOrq8D7g2yf30ruF3A7/W8VxrgU8kOdhRcuUKxCdJR7VU1aRjkCStoCQLwIauQ7m1c1zFeIack6Sp4m0VkjR7vgNcsJxJQIBfBBzrWNLcsedYkiRJauw5liRJkhqLY0mSJKmxOJYkSZIai2NJkiSp+T9/AK+PyGMNPgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Make a figure showing two signals with a beacon per signal\n",
    "fig, axes = plt.subplots(3,1, sharex=True, figsize=(12,4))\n",
    "axes[-1].set_xlabel(\"Time [ns]\")\n",
    "for i in range(0, 2):\n",
    "    axes[i].set_ylabel(\"Antenna {:d}\".format(i))\n",
    "    axes[i].plot(time/ns, impulses[i])\n",
    "    axes[i].plot(time/ns, beacons[i], marker='.')\n",
    "    axes[i].plot(time/ns, full_signals[i])\n",
    "\n",
    "# indicate timing of pulses\n",
    "for i, impulse_offset in enumerate(impulses_offsets):\n",
    "    kwargs = dict(color=['r','g'][i])\n",
    "    [ax.axvline(impulse_offset/ns, **kwargs) for ax in (axes[i], axes[-1])]\n",
    "\n",
    "# indicate timing of the beacons\n",
    "for i in range(0,2):\n",
    "    kwargs = dict(color=['r','g'][i], ls=(0, (3,1)))\n",
    "    tick_kwargs = dict(color='k', alpha=0.2)\n",
    "\n",
    "    tmp_beacon_offset = (beacon_init_phase / (2*np.pi*f_beacon) + i*beacon_time_offset) % (1/f_beacon)\n",
    "    # indicate every period of the beacon\n",
    "    beacon_ticks = [(n)*1/f_beacon + tmp_beacon_offset for n in range(1+int((time[-1] - time[0]) * f_beacon))]\n",
    "    \n",
    "    [axes[i].axvline(tick/ns, **{**kwargs, **tick_kwargs}) for tick in beacon_ticks]\n",
    "\n",
    "    # reference period in beacon\n",
    "    [ax.axvline(tmp_beacon_offset/ns, **kwargs) for ax in (axes[i], axes[-1])]\n",
    "    \n",
    "if not True:\n",
    "    axes[-1].set_xlim(0, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Solve it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "### correlation\n",
    "def correlation_and_lag(sig1, sig2, mode=\"full\", normalise=False):\n",
    "    corr = signal.correlate(sig1, sig2, mode=mode)\n",
    "    if normalise:\n",
    "        corr /= np.max(corr)\n",
    "\n",
    "    lags = signal.correlation_lags(sig1.size, sig2.size, mode=mode)\n",
    "    \n",
    "    return corr, lags\n",
    "\n",
    "def find_best_lag(sig1, sig2, fix_one_short=True, fix_positive=True, subtract_means=True, **corr_kwargs):\n",
    "    if subtract_means:\n",
    "        sig1 -= np.mean(sig1)\n",
    "        sig2 -= np.mean(sig2)\n",
    "    \n",
    "    corr, lags = correlation_and_lag(sig1, sig2, **corr_kwargs)\n",
    "    lag_id = corr.argmax()\n",
    "\n",
    "    lag = lags[lag_id]\n",
    "    \n",
    "    if fix_one_short:\n",
    "        # for some reason it is always one short\n",
    "        if lag > 0:\n",
    "            lag += 1\n",
    "        elif lag < 0:\n",
    "            lag -= 1\n",
    "\n",
    "    # turn negative lag into positive\n",
    "    if fix_positive and lag < 0:\n",
    "            lag += len(sig2)\n",
    "    \n",
    "    return lag, (corr, lags)\n",
    "\n",
    "def find_beacon_phase_delay(samplerate, f_beacon, reference_beacon, delayed_beacon, **lag_kwargs):\n",
    "    \"\"\"\n",
    "    Return phase delay of `beacon` with respect to `reference_beacon`.\n",
    "    Note that the returned value can be off by a multiple of $2\\pi$.\n",
    "    \n",
    "    Parameters\n",
    "    ==========\n",
    "    samplerate : float\n",
    "        Samplerate of both reference_beacon and delayed_beacon\n",
    "    f_beacon : float\n",
    "        Frequency of the beacons\n",
    "    reference_beacon : ndarray\n",
    "        The beacon to use as a reference\n",
    "    beacon : ndarray\n",
    "        The beacon to find the delay for\n",
    "    \"\"\"\n",
    "    \n",
    "    calc_lag, _ = find_best_lag(reference_beacon, delayed_beacon, **lag_kwargs)\n",
    "    \n",
    "    time = calc_lag / samplerate\n",
    "    \n",
    "    return 2*np.pi* f_beacon * time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# single out one period of the beacon\n",
    "beacon_samplerate = samplerate # MHz\n",
    "beacon_time = np.arange(0, 1/f_beacon, 1/beacon_samplerate)\n",
    "ref_beacon = sin_delay(f_beacon, beacon_time, phase=0, t_delay=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEWCAYAAABIVsEJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3xV9f3H8dcniySMQEhYCRBG2NsIKiAgQ0QER4uzgqPUKq5aV1u1WldrKy5c1dY9kGpFRdlgrYqEPUJI2GEl7Bkgyef3xznp7xqzSHJz7k0+z8fjPHLPvO9LQj75nvH9iqpijDHGnK4QrwMYY4wJTlZAjDHGVIgVEGOMMRViBcQYY0yFWAExxhhTIVZAjDHGVIgVEBPURORqEZlVRcdaICI3VmC/JBFREQmrihzBrjLfExH5o4i8U9WZjH9YATGecn9p7xeROhXZX1XfVdURVZ3LVJx9T2oPKyDGMyKSBAwEFBjjaRhTJawVVrtYATFeuhb4HngDGF/ahiIyQUQ2ishhEdkkIlf7LP/GZzsVkZtEJMNt2UwREXHXhYrI30Rkj3uMSaWdehKR60UkzT3OTBFpXcbnuV5EdojIThG5y+c4ISJyn4hsEJG9IjJVRGJ91n8kIrtE5KCIfC0iXX3WRbmZt7jrvxGRKHfdGBFZIyIH3JZcZ5/9NovIb0VkpbvfhyISWcq/7X9F5Hl323UiMtRnfYyIvO5+ru0i8qiIhBbZd7KI7AP+WMz35BwRWewee7GInOOzro2ILHS/r7OBuDL+jU0AsQJivHQt8K47nS8iTYvbSETqAs8BF6hqfeAcYHkpxx0NnAn0BMYB57vLfwlcAPQC+gAXl3QAEbkY+B1wKRAP/Ad4v4zPMwRIBkYA94nIMHf5be57DQJaAPuBKT77fenu1wRYivPvUeivwBk4nzkWuAcoEJEObp473HwzgM9EJMJn33HASKAN0AOYUEr2fsBGnF/gDwEf+xS5N4E8oD3Q2/18NxazbxPgMd+Dusf4Auf71xh4GvhCRBq7m7wHLHHf90+U8YeECTCqapNN1T4BA4BTQJw7vw64s4Rt6wIHgMuAqCLrJgDf+MwrMMBnfipwn/t6HvArn3XD3O3D3PkFwI3u6y+BG3y2DQGOAa2LyZfkHqeTz7K/AK+7r9OAoT7rmrufPayYYzV0jxXjvudxoGcx2z0ATC2Sbzsw2J3fDFxTJM/LJfz7TgB2AOKz7AfgF0BT4ITvvztwJTDfZ9+tJX1P3GP8UGT9d+42rXAKU12fde8B73j982lT+SZrgRivjAdmqeoed/49SvjrU1WPApcDNwE7ReQLEelUyrF3+bw+BtRzX7cAtvms831dVGvgWff00AFgHyBAQin7+B5vi/t+hcf6xOdYaUA+0NQ9rfake3rrEM4vfnD+Io8DIoENxbxXC/c9AFDVAvf9ffOV9O9QnO3q/gYvkr81EI7z716Y/xWc1kZxn7vUnD7HTnDX7Xe/v77rTJCwAmKqnXsOfxwwyD33vwu4E+gpIj2L20dVZ6rqcJy/3tcBf6/AW+8EEn3mW5ay7Tac1kpDnylKVb8tZR/f47XC+au+8FgXFDlWpKpuB64CxuK0hmJwWjPgFKs9QC7Qrpj32oHzy93Z2LnO0xKnFVIRCYXXiork34bTAonzyd5AVbv6bFtal94/yulz7O04349G7ilK33UmSFgBMV64GOcv8C441yN6AZ1xrjNcW3RjEWnqXjCui/PL7Ii7/+maCtwuIgki0hC4t5RtXwbuL7yg7V5I/nkZx39ARKLdfa4DPvQ51mOFF+FFJF5Exrrr6rufaS8QDTxeeDC3VfEP4GkRaeG2Vs4W55bnqcCFIjJURMKBu9zjlFbgStMEuE1Ewt3P2RmYoao7gVnA30SkgXtDQDsRGVTO484AOojIVSISJiKX43zfP1fVLUAq8LCIRIjIAOCiCuY3HrACYrwwHvinqm5V1V2FE/ACcHUxd0WF4PyC3IFzKmkQcHMF3vfvOL8MVwLLcH655VFMMVLVT4A/Ax+4p5ZW41yAL81CIBOYC/xVVQsfpnsWmA7MEpHDOHee9XPXvYVz2mY7sNZd5+u3wCpgMc5n/zMQoqrpwDXA8zgtlYuAi1T1ZJn/CsVbhHMhfw/OhfCfqeped921QISbbz8wDaclWCb3GKNxvn97cW4CGO1z6vIqnH+LfTgX79+qYH7jAfnxaU9jag8RuQDnwnJZt+fWaCIyAefmgQFeZzHBxVogptZwn6kY5Z5KScD5i/cTr3MZE6ysgJjaRICHcU7DLMO5G+pBTxMZE8TsFJYxxpgKsRaIMcaYCqlVHZ/FxcVpUlKS1zGMMSaoLFmyZI+qxhddXqsKSFJSEqmpqV7HMMaYoCIixfYQYKewjDHGVIgVEGOMMRViBcQYY0yFWAExxhhTIVZAjDHGVIinBURE/iEi2SKyuoT1IiLPiUimOzRnH59148UZtjRDRGwUM2OMqWZet0DewBlysyQX4PQQmgxMBF6C/w2T+RBOL559gYdEpJFfkxpjjPkRT58DUdWvRSSplE3GAm+5I6V9LyINRaQ5MBiYrar7AERkNk4hKmvM6gr5eGkWOw/mEl+vDnH1I4irV4fmMVHE16/jj7czxphKKyhQsvYfJyP7MOt3H2HCOUlERYRW6XsE+oOECfx4uMwsd1lJy39CRCbitF5o1apig519vnIn89Zl/2R5u/i6DGgfR//2cZzVrjENIsMrdHxjjKksVWXZtgN8vmInizbtZUPOEXJPFfxv/bkd4ujaIqZK3zPQC4gUs0xLWf7ThaqvAq8CpKSkVKjnyH9MOJPjJ/PZc+SEO51k054j/DdzL1NTs3jzuy2Ehggjuzbj+gFtOKO1nU0zxlSP9F2H+WTZdj5fuYOs/ceJCA2hb5tYru7Xmg5N69G+SX2Sm9bzyx+4gV5AsvjxONOJOKPSZeGcxvJdvsCfQaIiQmkZG03L2Gh3SVMmntuOE3n5LNt6gDlrd/Nh6ja+WLWT3q0acsOANozs2oywUK8vMxljaqKM3YeZPGc9M1btIixEGJAcxx3DOjCia9NqOxvieXfu7jWQz1W1WzHrLgQmAaNwLpg/p6p93YvoS4DCu7KWAmcUXhMpSUpKivqzL6yjJ/L4KHUb//x2M1v2HiO5ST0evbgb/do29tt7GmNqly17j/LMnAz+vXw70eGh3DCgDRP6tyG2boTf3lNElqhqyk+We1lAROR9nJZEHLAb586qcABVfVlEBGec7JHAMeA6VU11970e+J17qMdU9Z9lvZ+/C0ih/AJl5ppdPPZFGtsPHOdnZyRy/wWdaFzPLrobYyrmVH4BL87fwPPzMggLFcafncSvBrXza+EoFJAFpLpVVwEpdOxkHs/Py+TvX2+kXmQYvxvVmZ+fkYhTF40xpnw25BzhNx8uZ0XWQcb0bMEfLuxMkwaR1fb+JRWQQL8GEtSiI8K4d2QnLumdwB8+Wc0901ayaOM+HrukG5HhVXs7nTGm5ikoUN78bjNPfrmOqIhQXriqN6N7tPA61v9YAakGHZrW54OJZ/HcvAyenZtB2s5DvHzNGbRqHF32zsaYWunYyTzu+GA5s9bu5rxOTXjy0u7V2uooD7tFqJqEhAh3DOvAP8afSdb+Y4x+/j/ML+bZEmOM2X0ol3GvfMectN08MLoLr49PCbjiAVZAqt2QTk34/NaBJDaK5vo3F/PWd5u9jmSMCSBrdhxk7Av/ZVPOUV4bn8INA9oE7HVTKyAeaNU4mo9vPodhnZvy4KdrmDI/0+tIxpgAMG/dbn7+8neIwEc3ncN5nZp6HalUVkA8EhkeyotX92FsrxY8NTOdJ79cR226I84Y82Oz1uxi4ltLaBtfl3/f0p8uLRp4HalMdhHdQ+GhIUwe14t6dcJ4eeEGjpw4xSNjuhESEpjNVWOMf8xZu5tb3ltK14QY3r6hb9D0q2cFxGMhIcKjF3ejXmQYryzcSH6B8vgl3QP2nKcxpmrNW7ebX7+7hM7NG/DW9cFTPMAKSEAQEe4b2YlQEV5csIHGdevw2/M7eh3LGONn89OzuentpXRsVp+3r+9HTFTwFA+wAhIwRIS7z+/IvqMneWF+JnH1IpjQv43XsYwxfpK6eR+/ensJyU3r8c4N/YiJDq7iAVZAAoqIczpr39GTPPz5WmLr1WFMz8B56tQYUzU27znKL99KpUVMJG/f0I+G0f7vz8of7C6sABMWGsJzV/bmzKRY7pq6nK/X53gdyRhThfYfPcl1bywG4J/X9a2WzhD9xQpIAIoMD+Xv16bQLr4et7y7lMzsI15HMsZUgRN5+fzq7SVs33+cV69NoU1cXa8jVYoVkAAVExXO6xPOJCIshIlvp3Io95TXkYwxlaCq3DttJT9s3sdTP+/BmUmxXkeqNCsgASyhYRQvXt2HrXuPcecHyykosAcNjQlWLy3cwL+X7+C3IzowtleC13GqhBWQANevbWMevKgLc9dl88yc9V7HMcZUwHcb9vLXmelc2KM5twxp73WcKuNpARGRkSKSLiKZInJfMesni8hyd1ovIgd81uX7rJtevcmr1y/Oas24lESem5fJV6t3eh3HGHMasg/lcuv7y0iKq8ufL+tRox4S9uw2XhEJBaYAw4EsYLGITFfVtYXbqOqdPtvfCvT2OcRxVe1VXXm9JCI8MrYb63cf4TdTV9ChaX3axtfzOpYxpgx5+QVMem8ZR0/k8d4v+1GvTs16csLLFkhfIFNVN6rqSeADYGwp218JvF8tyQJQZHgoL13Th4iwEG59fxkn8vK9jmSMKcNTM9P5YfM+nri0Ox2a1vc6TpXzsoAkANt85rPcZT8hIq2BNsA8n8WRIpIqIt+LyMUlvYmITHS3S83JCe5nKprHRPGXy3qwZschnvoq3es4xphSzFqzi1e+3sg1Z7Xi4t4146J5UV4WkOJOBJZ0m9EVwDRV9f2zu5U7yPtVwDMi0q64HVX1VVVNUdWU+Pj4yiUOACO6NuMXZ7XmtW82sSDdRjQ0JhBlH8rlnn+tpFtCAx4Y3cXrOH7jZQHJAlr6zCcCO0rY9gqKnL5S1R3u143AAn58faRG+/2FnenYtD6//WgF2YdzvY5jjPGhqtw9bSW5p/J59ore1AkL9TqS33hZQBYDySLSRkQicIrET+6mEpGOQCPgO59ljUSkjvs6DugPrC26b00VGR7K81f15nBuHndNXWHPhxgTQN75fgsL1+fw+1GdaVfDb3bxrICoah4wCZgJpAFTVXWNiDwiImN8Nr0S+EB/PFxfZyBVRFYA84Enfe/eqg06NK3PA6O78J+MPbzx7Wav4xhjgA05R3hsRhqDOsRzzVmtvY7jd1KbhlFNSUnR1NRUr2NUGVXlhjdT+XbDHr66/VySgrxfHWOC2an8Ai576Vu27jvGrDvOpUmDSK8jVRkRWeJec/4RexI9iIkIj1/SnfDQEO6ZttJOZRnjoefnZbIy6yBPXNK9RhWP0lgBCXLNYiJ5cHQXfti8j7e+2+x1HGNqpdXbDzJlfiaX9knggu7NvY5TbayA1AA/OyORIR3j+fNX6WzZe9TrOMbUKnn5Bdz38Upi60bw0OiuXsepVlZAagAR4YlLexAWKtxtp7KMqVavfbOJ1dsP8ciYrkE5LG1lWAGpIZrFRPLA6C78sGkfb3+/xes4xtQKm/ccZfLs9Yzo0pSR3Zp5HafaWQGpQX5+RiLndojnqZnp7DpoDxga40+qyv0fryIiLIQ/XdytRvWyW15WQGoQEeHRsd04lV/AI5+v8TqOMTXah4u38d3GvfxuVGea1pK7roqyAlLDtGoczW1Dk5mxahfz11lfWcb4Q/bhXB6bkUa/NrFcntKy7B1qKCsgNdAvB7alfZN6PDh9NcdPWrfvxlS1J2as48SpAp64tDshIbXv1FUhKyA1UERYCI9e3I1t+47zwvwMr+MYU6Ms2riXT5Zt51eD2tb6gd2sgNRQZ7VtzGV9Enn1641k7D7sdRxjaoRT+QU8+OkaEhpGcfPgmjO2eUVZAanBfjeqE9ERYfz+36upTX2eGeMvb367mfTdh3nooi5ERdTcbtrLywpIDda4Xh3uHdmJHzbt47OVO72OY0xQ230ol2fmZDCkYzzDuzT1Ok5AsAJSw11+Zku6tmjAEzPSOHYyz+s4xgStx2ekcTK/gD+O6Vorn/kojhWQGi40RPjjmK7sPJjLyws2eB3HmKD03Ya9fLp8BzcNakfrxjZsQiErILXAmUmxjOnZgle+3si2fce8jmNMUMkvUB7+rPDCeTuv4wQUTwuIiIwUkXQRyRSR+4pZP0FEckRkuTvd6LNuvIhkuNP46k0efO4f1YkQER6fkeZ1FGOCytTUbazbdZjfjepMZLhdOPflWQERkVBgCnAB0AW4UkS6FLPph6ray51ec/eNBR4C+gF9gYdEpFE1RQ9KzWOcv56+XL2Lbzfs8TqOMUHhcO4p/jYrnTOTGjGqe+3rLLEsXrZA+gKZqrpRVU8CHwBjy7nv+cBsVd2nqvuB2cBIP+WsMX55blsSG0Xx8PS15OUXeB3HmIA3Zf4G9hw5yQOju9iF82J4WUASgG0+81nusqIuE5GVIjJNRAo7nSnvvojIRBFJFZHUnJycqsgdtCLDQ/nDhZ1J332YD1O3lb2DMbXY1r3H+Mc3m7i0TwI9Eht6HScgeVlAiivnRZ92+wxIUtUewBzgzdPY11mo+qqqpqhqSnx8fIXD1hTnd21G36RYJs/O4MgJu63XmJI8+VUaoSHCPed38jpKwPKygGQBvt1YJgI7fDdQ1b2qesKd/TtwRnn3NcUTEX53YWf2HDnBqwvttl5jirNo415mrNrFTYPa0SymdnbVXh5eFpDFQLKItBGRCOAKYLrvBiLiOzr9GKDwFqKZwAgRaeRePB/hLjPl0KtlQy7q2YJX/7PRBp4ypoiCAuXRL9JoHhPJxHPbeh0noHlWQFQ1D5iE84s/DZiqqmtE5BERGeNudpuIrBGRFcBtwAR3333An3CK0GLgEXeZKad7zu9IQQE8PTvd6yjGBJTPVu5g1faD/HZER+vvqgxSmzrZS0lJ0dTUVK9jBIzHvljLa99sYsZtA+ncvIHXcYzx3Im8fIY9vZC6EWF8cdtAQmvxWB++RGSJqqYUXW5Potdik4Yk0yAy3B4uNMb17vdb2bbvOPeP6mzFoxysgNRiMdHh3DY0mf9k7GHh+tp9i7Mxh3JP8fy8DPq3b8y5yXFexwkKVkBquV+c1ZpWsdH8+ct1FBTUntOZxhT1ysIN7D92ivtGdraHBsvJCkgtFxEWwl0jOrB25yE+X2VjhpjaadfBXF7/ZhNjerage2KM13GChhUQw0U9WtCpWX3+NiudU9bFiamFJs9eT36Bcvf5Hb2OElSsgBhCQoR7R3Ziy95jfLDYujgxtUtm9mE+WrKNa85qTcvYaK/jBBUrIAaAwR3j6ZsUy3NzM2zkQlOrPD17PVHhoUwa0t7rKEHHCogBnC5O7hnZkZzDJ/jnfzd7HceYarEq6yAzVu3ihoFtaVyvjtdxgo4VEPM/KUmxDOvchJcXbuDAsZNexzHG7/46K52G0eHcOLCN11GCkhUQ8yO/Pb8jR07k8ZKNn25quEUb97JwfQ6/HtSOBpHhXscJSlZAzI90ataAS3ol8Ma3m8k+ZB0tmppJVfnrrHSaNqjD+HOSvI4TtKyAmJ+4fVgy+QXKi9YKMTXUgvQcFm/ez63nJds455VgBcT8ROvGdfl5SkveW7SV7QeOex3HmCpVUKA8NTOdVrHRjEtpWfYOpkRWQEyxbj3PuaXxhXkZHicxpmp9uXoXa3ce4s7hyUSE2a/AyrB/PVOsFg2juKpfK6amZrF5z1Gv4xhTJfILlMlz1tO+ST3G9EzwOk7Q87SAiMhIEUkXkUwRua+Y9b8RkbUislJE5opIa591+SKy3J2mF93XVN7Ng9sRHio8N9daIaZm+GzFDjKzj3DnsA7WXXsV8KyAiEgoMAW4AOgCXCkiXYpstgxIUdUewDTgLz7rjqtqL3cag6lyTRpEMv7sJD5Zvp2M3Ye9jmNMpeTlF/Ds3Aw6NavPBd2aeR2nRvCyBdIXyFTVjap6EvgAGOu7garOV9Vj7uz3QGI1Z6z1fjWoHdHhoTwzx1ohJrh9smw7m/Yc5c7hHQix1keV8LKAJAC+PfdluctKcgPwpc98pIikisj3InJxSTuJyER3u9ScHBs06XTF1o3ghgFt+GLVTtbuOOR1HGMq5FR+Ac/Ny6BbQgNGdGnqdZwaw8sCUtyfAMWOaCQi1wApwFM+i1u5Y/ReBTwjIu2K21dVX1XVFFVNiY+Pr2zmWumGgW2pHxnGs3PXex3FmAqZtiSLbfuO85vhHWywqCrkZQHJAnxvwk4EdhTdSESGAb8HxqjqicLlqrrD/boRWAD09mfY2iwmKpwbBrRh5prdrN5+0Os4xpyWE3n5PD83g14tGzKkYxOv49QoZRYQEYkWkQdE5O/ufLKIjK6C914MJItIGxGJAK4AfnQ3lYj0Bl7BKR7ZPssbiUgd93Uc0B9YWwWZTAmuH9CGBpFhdi3EBJ2pi7ex42Aud42w1kdVK08L5J/ACeBsdz4LeLSyb6yqecAkYCaQBkxV1TUi8oiIFN5V9RRQD/ioyO26nYFUEVkBzAeeVFUrIH7UIDKcGwe2ZU7ablZlWSvEBIfcU/lMmb+BM5MaMaB9nNdxapywcmzTTlUvF5ErAVT1uFRRGVfVGcCMIsse9Hk9rIT9vgW6V0UGU37X9U/i9W828cyc9bw+4Uyv4xhTpqmp29h1KJenx/W01ocflKcFclJEonAvcLsXq0+UvoupiepHhvPLgW2Yuy6bFdsOeB3HmFLlnsrnxfkb6JsUy9ntGnsdp0YqTwF5CPgKaCki7wJzgXv8msoErPHnJNEwOpxn5tgdWSawFbY+7hiWbK0PPymzgKjqbOBSYALwPs6T4Qv8G8sEKqcV0pb56Tks27rf6zjGFMu59pFprQ8/K7GAiEifwgloDezEuc22lbvM1FLjz0miUXQ4z1ofWSZAfbh4G7sPnbDWh5+VdhH9b+7XSJyH+FbgPPzXA1gEDPBvNBOo6tUJ48aBbXlqZjrLtx2gV8uGXkcy5n9yT+Xz4gJrfVSHElsgqjpEVYcAW4A+7tPcZ+A8sJdZXQFNYCq8FmI99ZpAY62P6lOei+idVHVV4YyqrgZ6+S+SCQb16oRx44A2zFuXzcosuyPLBAZrfVSv8hSQNBF5TUQGi8gg94n0NH8HM4Fv/DlJxERZK8QEjqmpTuvjdmt9VIvyFJDrgDXA7cAdOF2GXOfPUCY41I90+siak5ZtfWQZz53Iy+elBRtIad2Ic6z1US3KcxtvrqpOVtVL3GmyquZWRzgT+Cb0T6JBZJjdkWU891FqFjsP5lrroxqVpzPFTSKysehUHeFM4GsQGc4NA9oye+1u1uywVojxxsm8Al5asIHerRpan1fVqDynsFKAM91pIPAc8I4/Q5ngMqF/EvUjw+xaiPHMv5Zmsf3AcW4faq2P6lSeU1h7fabtqvoMcF41ZDNBIiYqnOv7O+OFpO20UQtN9TqVX8CU+Zn0bNmQQR1s0LjqVJ5TWH18phQRuQmoXw3ZTBC5vn8b6tcJ4/l51gox1euTpdvJ2n+c24e2t9ZHNStPd+5/83mdB2wCxvknjglWMdHhTOifxPPzMknfdZiOzexvDON/p/ILeGF+Jt0TYmy0QQ+U5xrIDYVPpavqcFWdCJz0dzATfK7v34a6EaG8MN86KjDV49PlO9i67xi32bUPT5SngEwr57LTJiIjRSRdRDJF5L5i1tcRkQ/d9YtEJMln3f3u8nQROb8q8pjKaVQ3gvHnJPH5yh1kZh/2Oo6p4fLcax9dmjdgWGdrfXihtN54O4nIZUCMiFzqM03A6WCxUkQkFJgCXAB0Aa4UkS5FNrsB2K+q7YHJwJ/dfbvgjKHeFRgJvOgez3jsxoFtiQoP5YV51gox/vXZyh1s2nPUWh8eKq0F0hEYDTQELvKZ+gC/rIL37gtkqupGVT0JfACMLbLNWOBN9/U0YKg7nO5Y4ANVPaGqm3A6d+xbBZlMJcXWjeAXZ7Vm+oodbMw54nUcU0PlFyjPz8ukU7P6jOjS1Os4tVZpvfF+qqrXAaNV9Tqf6TZ3TPLKSgC2+cxnucuK3UZV84CDQONy7guAiEwUkVQRSc3JyamC2KYsvzy3LRFhIXYtxPjNF6t2sjHHaX2EhFjrwyulncIqHLb2KhF5ruhUBe9d3Hddy7lNefZ1Fqq+6nZFnxIfb/eIV4e4enW4pl9r/r1sO5v3HPU6jqlhCgqU5+dmkNykHiO7NvM6Tq1W2imswh53U4ElxUyVlQW09JlPxBnxsNhtRCQMiAH2lXNf46GJg9oSHhrCFGuFmCr25epdZGQf4VZrfXiuxOdAVPUz9+ubJW1TSYuBZBFpA2zHuSh+VZFtpgPjge+AnwHzVFVFZDrwnog8DbQAkoEf/JTTVECT+pFc1a8Vb323hduGJtMyNtrrSKYGKChQnp+XQdv4ulzYvbnXcWq9EguIiHxGCaeFAFR1TGXeWFXzRGQSMBMIBf6hqmtE5BEgVVWnA68Db4tIJk7L4wp33zUiMhWna/k84BZVza9MHlP1bhrUjncXbWXK/EyevKyH13FMDTBr7S7W7TrM5Mt7EmqtD8+JavE1QkQGlbajqi70SyI/SklJ0dTUVK9j1CoPfrqa9xZtZcHdg0lsZK0QU3GqyqjnvuH4yTzm/GYQYaHleYzNVAURWaKqKUWXl3YX1sLCCecU0n6cVsB3wVg8jDd+PbgdISK8tGCD11FMkJu91umsc9J5yVY8AkR5OlO8ENiA0437C0CmiFzg72CmZmgeE8XPUxKZmrqNHQeOex3HBClV5dm5GbRuHM3FvVp4Hce4ylPG/wYMUdXBqjoIGILzVLgx5fLrwe0AeHmhtUJMxcxNy2bNjkPcMqS9tT4CSHm+E9mq6nsv5kYg2095TA2U2Cian52RyAc/bGPXQRsN2ZyewtZHy9goLuld7PPCxiPlKSBrRGSGiEwQkfHAZ8Diwr6x/JzP1BA3D25PgeL6B5gAABqPSURBVKq1QsxpW5Cew6rtB5k0pD3h1voIKOX5bkQCu4FBwGAgB4jF6RdrtN+SmRqlZWw0l/VJ5L0ftrL7kLVCTPmoKs/MzSChYRSX9kn0Oo4poswBpdz+sIyptFuGtGfa0ixeWrCBP47p6nUcEwQWrs9hxbYDPH5Jd2t9BKAyC4j7pPitQJLv9pV9kNDUPq0aR3NZnwTe+2Ervx7cjqYNKj0qgKnBCq99JDSM4mdnWOsjEJWnpP8b2Aw8j3NHVuFkzGmbNCSZ/AK7FmLK9nXGHpZtPcCvB7cjIsxaH4GoPGOi56pqVfS+awytGkdzae8E3lu0lV8PakcTa4WYYqgqk2evJ6FhFONSWpa9g/FEecr6syLykIicLSJ9Cie/JzM11qTz2pNXoLy8cKPXUUyAWrg+h+XbDnDLkPbW+ghg5WmBdAd+AZwHFLjL1J035rS1blyXS3on8O6iLdw0qK21QsyPqCrPzLFrH8GgPKX9EqCtqg5S1SHuZMXDVMqkIU4r5CW7FmKKWGCtj6BRnu/OCpxx0Y2pMklxdbm0dwLvLrLnQsz/s9ZHcClPAWkKrBORmSIy3Z0+9XcwU/Pdel4yBQXKizZqoXEtSHee+5h0nrU+gkF5roE85PNagAHAlf6JY2qTVo2dPrLe/2EbNw1uR/OYKK8jGQ85rQ/nzqvL7KnzoFBmiXfH/jgIXAi8AQwFXq7Mm4pIrIjMFpEM92ujYrbpJSLficgaEVkpIpf7rHtDRDaJyHJ36lWZPMY7twxpj6I2drph3rpsVmQdtNZHECnxuyQiHUTkQRFJwxkHZBvOCIZDVPX5Sr7vfcBcVU0G5rrzRR0DrlXVrsBI4BkR8b0Wc7eq9nKn5ZXMYzzSMjaacSkt+XDxNrL2H/M6jvGIqvL07PW0io22ax9BpLQyvw6ntXGRqg5wi0ZVjTs+FnjTff0mcHHRDVR1vapmuK934HQhH19F728CyC1D2iOItUJqsZlrdrFmxyFuG5psfV4FkdK+U5cBu4D5IvJ3ERmKcw2kKjRV1Z0A7tcmpW0sIn2BCJyREQs95p7amiwidUrZd6KIpIpIak5OTlVkN1WsRcMorujbko9Ss9i2z1ohtU1BgTJ5dgZt4+raaINBprQx0T9R1cuBTsAC4E6gqYi8JCIjyjqwiMwRkdXFTGNPJ6CINAfeBq5T1cIHGe93c52J07X8vaV8jldVNUVVU+LjrQETqG4e3J6QEOG5uRleRzHV7PNVO0nffZjbh9lY58GmPBfRj6rqu6o6GkgEllP8NYui+w1T1W7FTJ8Cu93CUFggih3hUEQaAF8Af1DV732OvVMdJ4B/An3L8VlNAGsWE8nV/Vrx8bLtbMw54nUcU03y8gt4Zs56OjStx0U9rPURbE6r3KvqPlV9pQqeRJ8OjHdfjwd+8lyJiEQAnwBvqepHRdYVFh/BuX6yupJ5TAC4eXB7IkJDmDzHWiG1xafLd7Ax5yh3DutASEhVnSE31cWr9uKTwHARyQCGu/OISIqIvOZuMw44F5hQzO2674rIKmAVEAc8Wr3xjT/E16/Ddf2T+GzFDtJ2HvI6jvGzU/kFPDs3gy7NG3B+12ZexzEVUJ4HCaucqu7FucOr6PJU4Eb39TvAOyXsb31x1VC/Orcdb3+/hb/NWs9r41O8jmP86F9Lsti67xivXZtirY8gZVesTECJiQ5n4sC2zEnbzbKt+72OY/wk91Q+z87NoHerhgztXOpNmCaAWQExAee6AW2IrRvB32at9zqK8ZN3vt/CzoO53H1+R5xLmSYYWQExAadenTBuHtyObzL38N2GvV7HMVXscO4ppszPZGByHOe0i/M6jqkEKyAmIF1zVmuaNqjDX2elo6pexzFV6PVvNrH/2CnuPr+j11FMJVkBMQEpMjyU24Yms2TLfuatK/YxIROE9h09yWv/2cQF3ZrRI9GGGQp2VkBMwBqX0pI2cXX5y1fp5BdYK6QmeHF+JsdO5nHXiA5eRzFVwAqICVjhoSHcNaID6bsP88my7V7HMZW048Bx3vp+C5f2SaR9k/pexzFVwAqICWijujWnR2IMk2evJ/dUVXUGbbzw3NwMULhjWLLXUUwVsQJiAlpIiHDvyE5sP3Ccd77f4nUcU0EZuw8zNXUbV5/VisRG0V7HMVXECogJeP3bxzEwOY4X5mdyKPeU13FMBfz5q3XUrRPGbedZ66MmsQJigsK9Iztx4NgpXlm4oeyNTUD5fuNe5qRlc/Pg9jSqG+F1HFOFrICYoNAtIYYxPVvw+jeb2H0o1+s4ppxUlSdmpNE8JpLr+id5HcdUMSsgJmjcNaID+QXK5NnWxUmw+GLVTlZkHeSuER2JDA/1Oo6pYlZATNBo3bgu156dxIep21i7w7p7D3Qn8wr4y1fpdGpWn0t6J3gdx/iBFRATVG47L5mYqHAen5FmXZwEuHcXbWHrvmPcP6ozodZde43kSQERkVgRmS0iGe7XRiVsl+8zmNR0n+VtRGSRu/+H7uiFphaIiQ7ntvOS+SZzDwvSc7yOY0pw8NgpnpubQf/2jTk32TpMrKm8aoHcB8xV1WRgLiWPsX5cVXu50xif5X8GJrv77wdu8G9cE0iuOas1beLq8ugXazmVX+B1HFOMZ+dmcOD4KX4/qot1116DeVVAxgJvuq/fxBnXvFzccdDPA6ZVZH8T/CLCQrj/gk5syDnKBz9s9TqOKSIz+whvfbeZK85sSZcWDbyOY/zIqwLSVFV3ArhfSxqSLFJEUkXkexEpLBKNgQOqmufOZwF2ha6WGd6lKf3axDJ5ToY9XBhgHp+RRlR4KHeNsO7aazq/FRARmSMiq4uZxp7GYVqpagpwFfCMiLQDimsPl3g1VUQmukUoNSfHzpnXFCLCA6O7sP/YSabMy/Q6jnEtXJ/DvHXZ3Dq0PXH16ngdx/iZ3wqIqg5T1W7FTJ8Cu0WkOYD7tdgBH1R1h/t1I7AA6A3sARqKSJi7WSKwo5Qcr6pqiqqmxMfHV9nnM97rlhDDZX0S+cd/N7Eh54jXcWq9vPwC/vT5WpIaRzPhnDZexzHVwKtTWNOB8e7r8cCnRTcQkUYiUsd9HQf0B9aqc+/mfOBnpe1vaod7R3YiMiyUhz9ba7f1euzdRVvJzD7C70Z1JiLMnhCoDbz6Lj8JDBeRDGC4O4+IpIjIa+42nYFUEVmBUzCeVNW17rp7gd+ISCbONZHXqzW9CRjx9etwx/AOfL0+h9lrd3sdp9Y6cOwkk+esp3/7xgzv0tTrOKaahJW9SdVT1b3A0GKWpwI3uq+/BbqXsP9GoK8/M5rgce3Zrfngh6386Yu1nNsh3rrM8MBfZqZzODePB0bbbbu1ibUzTdALDw3h4TFd2bbvOK9+vdHrOLXO8m0HeP+HrUw4J4lOzey23drECoipEc5pH8eF3ZszZX4mWfuPeR2n1sgvUP7w71U0qV/HRhqshayAmBrjdxd2RgQe/TzN6yi1xnuLtrB6+yH+cGEX6keGex3HVDMrIKbGSGgYxaQh7flqzS7mrbML6v6Wc/gEf5mZzoD2cYzu0dzrOMYDVkBMjTLx3HYkN6nHA/9ew9ETeWXvYCrsiS/TyD2Vz8Nju9qF81rKCoipUSLCQnj80u5sP3DcBp7yo0Ub9/Lx0u1MPLct7eLreR3HeMQKiKlxzkyK5ap+rfjHfzexevtBr+PUOLmn8rn/k1UkNopi0hC7cF6bWQExNdK9IzvRuF4d7vt4JXnW5XuVem5uBhtzjvLEpd2JirBnbmozKyCmRoqJCuehi7qwevsh3vh2s9dxaozV2w/yytcbGZeSyMBk61uutrMCYmqsC7s3Z0jHeJ6evZ5t++zZkMo6lV/APdNWEls3gt+P6uJ1HBMArICYGktE+NPF3RDgnmkrKSiwzhYr49WvN7J25yH+NLYbMdH2zIexAmJquMRG0TwwugvfbdzLW99t9jpO0MrMPsKzczK4sHtzRnZr5nUcEyCsgJga7/IzWzKkYzxPfrWOjTZuyGnLyy/gnmkriK4Tyh/HdPU6jgkgVkBMjSciPHlZD+qEhXLXRyvsrqzTNGX+BpZuPcDDY7oSX99GGTT/zwqIqRWaNojkkbFdWbb1AK/+x3rsLa+lW/fz3LwMLu7VgrG9EryOYwKMFRBTa4zp2YJR3ZsxefZ60nYe8jpOwDtyIo87P1xOswaRPHJxN6/jmADkSQERkVgRmS0iGe7XRsVsM0RElvtMuSJysbvuDRHZ5LOuV/V/ChNsRIRHL+5OTFQEt76/jGMnra+s0jzy2Rq27jvG0+N60sB62jXF8KoFch8wV1WTgbnu/I+o6nxV7aWqvYDzgGPALJ9N7i5cr6rLqyW1CXqxdSN45vJebMg5woOfrvE6TsD6avVOpqZm8etB7ejXtrHXcUyA8qqAjAXedF+/CVxcxvY/A75UVXsazFTagOQ4Jg1pz7QlWfxrSZbXcQLO9gPHue/jVXRPiOGOYR28jmMCmFcFpKmq7gRwvzYpY/srgPeLLHtMRFaKyGQRKfHWEBGZKCKpIpKak5NTudSmxrh9aDJ928TywKerycy2W3sLncjL5+Z3lpCXrzx7RS8iwuwyqSmZ3346RGSOiKwuZhp7msdpDnQHZvosvh/oBJwJxAL3lrS/qr6qqimqmhIfb333GEdYaAjPXdGbyPBQJr23lNxT+V5HCgiPfLaWFVkH+evPe9LWumk3ZfBbAVHVYararZjpU2C3WxgKC0R2KYcaB3yiqqd8jr1THSeAfwJ9/fU5TM3VLCaSp8f1ZN2uwzz06RpUa3dXJ9OWZPHuoq3cNKidPW1uysWr9ul0YLz7ejzwaSnbXkmR01c+xUdwrp+s9kNGUwsM7tiESUPa82HqNt6sxb32rtlxkN9/soqz2zbmtyPsuocpH68KyJPAcBHJAIa784hIioi8VriRiCQBLYGFRfZ/V0RWAauAOODRashsaqjfDO/A8C5NeeTztXy9vvZdJzt47BS/fmcpjaIjeP6q3oSF2nUPUz5Sm5rtKSkpmpqa6nUME4COnsjjspe+ZfuB43xyc3/aN6kd5/9P5OVz7es/sHTrfj781dn0afWTR7KMQUSWqGpK0eX2p4YxQN06Ybw2PoU6YSHc+OZiDhw76XUkvysoUO7+aCWLNu3jrz/vacXDnDYrIMa4EhtF88ovzmDHgVxufncpJ/NqdqeLT81KZ/qKHdwzsqP1c2UqxAqIMT7OaB3LE5d259sNe7njw2U1tufed77fwksLNnBVv1b8elA7r+OYIBXmdQBjAs1lZySy/9hJHv0ijeiIVfzlsh6EhIjXsarM7LW7efDT1ZzXqQmPjOmKczOjMafPCogxxbhxYFuOnMjjmTkZ1I1wBlKqCb9o56bt5pZ3l9ItIYbnr7Q7rkzlWAExpgS3D03m6Ik8/v6fTdStE8Y9Izt5HalSZq3ZxS3vLaVz8wa8fX0/6tax//6mcuwnyJgSiAi/G9WZoyfzeXHBBhS45/yOQdkS+Wr1Lia9t5SuCTG8dX1fYqKse3ZTeVZAjCmFiPDoWGcwpZcWbCDn8AmeuLQ74UF06ufLVTu59f1ldE+M4c3r+9rYHqbKWAExpgwhIcJjF3ejaf1IJs9Zz54jJ3jx6j5ERwT2fx9V5fVvNvH4jDR6t2rEG9edSX0rHqYKBc+fUcZ4SES4fVgyT1zana/X53Dl3xex98gJr2OV6ERePvdMW8mjX6RxftdmvH1DXysepspZATHmNFzZtxWv/CKFdTsPMeaF/5K6eZ/XkX5iz5ETXP33RXy0JIvbhiYz5arAby2Z4GQFxJjTNLxLUz666WzCQoXLX/2eKfMzyS8IjD7lFm/ex9gX/suq7Qd5/sre/GZ4hxr1DIsJLFZAjKmAHokN+fzWAYzq3pynZqbzi9cXkX0o17M8x07m8fBnaxj3yneIwEc3nc1FPVt4lsfUDtYbrzGVoKp8lJrFg9NXExkeyh1Dk7n6rNbVepfWoo17uedfK9my9xjXnt2ae0d2smc8TJUqqTdeKyDGVIHM7MM8+Okavt2wl3bxdfn9hZ0Z0rGJX58Z2ZBzhBfmZfLJsu20io3mz5f14Ox2jf32fqb2sgKCFRDjX6rKnLRsHp+RxqY9RxmYHMeNA9sysH1clV6HWL/7MC/My+SzlTuIDAvl2nNac/vQZLtQbvwmoAqIiPwc+CPQGeirqsX+VheRkcCzQCjwmqoWjlzYBvgAiAWWAr9Q1TIHcLACYqrDybwC3vl+Cy/Mz2Tf0ZMkNIxiXEpLxp2ZSPOYqAodc9/Rk8xJ281Xq3cxPz2bqPBQrj07iRsHtiGuXp0q/gTG/FigFZDOQAHwCvDb4gqIiIQC63GGvM0CFgNXqupaEZkKfKyqH4jIy8AKVX2prPe1AmKq04m8fGav3c0HP2zjm8w9hAh0S4ihV8uG9G7VkF4tG5HUOPonp7lO5Rewff9xNu89Smb2EeamZfPD5n3kFygtYiK5tE8i1w9oQ2zdCI8+maltAqqA/O/NRRZQcgE5G/ijqp7vzt/vrnoSyAGaqWpe0e1KYwXEeGXr3mP8a2kWizbtZWXWQY6dzAcgNESIDAshKiKUyPBQVGHXodwf3Rbcvkk9RnZtxvldm9EtoUFQ9sVlgltJBSSQT5omANt85rOAfkBj4ICq5vksL3E4NRGZCEwEaNWqlX+SGlOGVo2juXN4BwDyC5SM7MMs23qAbfuOkXuqgNy8fHJP5lOgSmKjaFo3jiYpri6tG0fTpH6kx+mNKZ7fCoiIzAGaFbPq96r6aXkOUcwyLWV5sVT1VeBVcFog5XhfY/wqNETo1KwBnZo18DqKMZXitwKiqsMqeYgsoKXPfCKwA9gDNBSRMLcVUrjcGGNMNQrkJ9EXA8ki0kZEIoArgOnqXLSZD/zM3W48UJ4WjTHGmCrkSQERkUtEJAs4G/hCRGa6y1uIyAwAt3UxCZgJpAFTVXWNe4h7gd+ISCbONZHXq/szGGNMbWcPEhpjjClVSXdhBfIpLGOMMQHMCogxxpgKsQJijDGmQqyAGGOMqZBadRFdRHKALRXcPQ7nGZRgYFn9w7JWvWDJCbU7a2tVjS+6sFYVkMoQkdTi7kIIRJbVPyxr1QuWnGBZi2OnsIwxxlSIFRBjjDEVYgWk/F71OsBpsKz+YVmrXrDkBMv6E3YNxBhjTIVYC8QYY0yFWAExxhhTIVZAykFERopIuohkish9XufxJSL/EJFsEVntsyxWRGaLSIb7tZGXGd1MLUVkvoikicgaEbk9gLNGisgPIrLCzfqwu7yNiCxys37oDjMQEEQkVESWicjn7nxAZhWRzSKySkSWi0iquyzgfgYARKShiEwTkXXuz+3ZgZhVRDq6/56F0yERuaM6sloBKYOIhAJTgAuALsCVItLF21Q/8gYwssiy+4C5qpoMzHXnvZYH3KWqnYGzgFvcf8dAzHoCOE9VewK9gJEichbwZ2Cym3U/cIOHGYu6HWfYg0KBnHWIqvbyeU4hEH8GAJ4FvlLVTkBPnH/fgMuqqunuv2cv4AzgGPAJ1ZFVVW0qZcIZs2Smz/z9wP1e5yqSMQlY7TOfDjR3XzcH0r3OWEzmT4HhgZ4ViAaWAv1wnuwNK+7nwuOMie4viPOAz3GGfQ7UrJuBuCLLAu5nAGgAbMK90SiQsxbJNwL4b3VltRZI2RKAbT7zWe6yQNZUVXcCuF+beJznR0QkCegNLCJAs7qnhJYD2cBsYANwQJ2BziCwfg6eAe4BCtz5xgRuVgVmicgSEZnoLgvEn4G2QA7wT/fU4GsiUpfAzOrrCuB997Xfs1oBKZsUs8zufa4gEakH/Au4Q1UPeZ2nJKqar84pgUSgL9C5uM2qN9VPichoIFtVl/guLmZTz7O6+qtqH5xTwreIyLleBypBGNAHeElVewNHCYDTVaVxr3ONAT6qrve0AlK2LKClz3wisMOjLOW1W0SaA7hfsz3OA4CIhOMUj3dV9WN3cUBmLaSqB4AFONdtGopImLsqUH4O+gNjRGQz8AHOaaxnCMysqOoO92s2znn6vgTmz0AWkKWqi9z5aTgFJRCzFroAWKqqu915v2e1AlK2xUCye1dLBE4TcbrHmcoyHRjvvh6Pc73BUyIiOGPXp6nq0z6rAjFrvIg0dF9HAcNwLqDOB37mbhYQWVX1flVNVNUknJ/Neap6NQGYVUTqikj9wtc45+tXE4A/A6q6C9gmIh3dRUOBtQRgVh9X8v+nr6A6snp90ScYJmAUsB7nPPjvvc5TJNv7wE7gFM5fTTfgnAOfC2S4X2MDIOcAnNMoK4Hl7jQqQLP2AJa5WVcDD7rL2wI/AJk4pwnqeJ21SO7BwOeBmtXNtMKd1hT+XwrEnwE3Vy8g1f05+DfQKICzRgN7gRifZX7Pal2ZGGOMqRA7hWWMMaZCrIAYY4ypECsgxhhjKsQKiDHGmAqxAmKMMaZCrIAYY4ypECsgxpSTiDT26TJ7l4hs95n/1g/vN0FEckTktQru/5Sb87dVnc0YcPp7McaUg6ruxXm4DBH5I3BEVf/q57f9UFUnVWRHVb1bRI5WdSBjClkLxJgqICJH3K+DRWShiEwVkfUi8qSIXO0OULVKRNq528WLyL9EZLE79S/He0wQkY9F5Ct3kKC/uMtDReQNEVntvsed/v20xjisBWJM1euJ03vvPmAj8Jqq9hVnFMZbgTtwBiuarKrfiEgrYCbF9/hbVC+crvBPAOki8jxON90JqtoNnJH0qvoDGVMcKyDGVL3F6o7DICIbgFnu8lXAEPf1MKCL08ckAA1EpL6qHi7j2HNV9aB77LVAa5x+pdq6xeQLn/czxq+sgBhT9U74vC7wmS/g///PhQBnq+rxShw7H2fUwf0i0hM4H7gFGAdcf9qpjTlNdg3EGG/MAv53cVxEelX0QCISB4So6r+AB3DGrTDG76wFYow3bgOmiMhKnP+HXwM3VfBYCThDrxb+QXh/FeQzpkzWnbsxAUpEJgApFb2N1z3GH6me241NLWSnsIwJXMeBCyrzICFwDc543sZUOWuBGGOMqRBrgRhjjKkQKyDGGGMqxAqIMcaYCrECYowxpkL+D0PP1rYNRqRGAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# show beacon period\n",
    "fig, ax = plt.subplots()\n",
    "ax.set_title(\"A single beacon period\")\n",
    "ax.set_xlabel(\"Time [ns]\")\n",
    "ax.set_ylabel(\"Amplitude\")\n",
    "ax.plot(ref_beacon)\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#####  1.1 Beacon Phase Delays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 0.0 0.0\n",
      "37 0.0074 3.2546899891190257\n",
      "Beacon delays [ns] \\pm k*14.285714285714285ns: [0.  7.4]\n"
     ]
    }
   ],
   "source": [
    "beacon_phase_delays = np.array([\n",
    "    find_beacon_phase_delay(beacon_samplerate, f_beacon, beacons[0], beacon)\n",
    "        for beacon in beacons\n",
    "])\n",
    "\n",
    "beacon_time_delays = (beacon_phase_delays) / (2*np.pi * f_beacon)\n",
    "\n",
    "print(\"Beacon delays [ns] \\pm k*{}ns: {}\".format(1/f_beacon/ns, beacon_time_delays/ns))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make a figure showing the corrected beacons\n",
    "fig, ax = plt.subplots(1,1, sharex=True)\n",
    "ax.set_xlabel(\"Time [ns]\")\n",
    "ax.set_ylabel(\"Amplitude [au]\")\n",
    "ax.set_title(\"Beacon delays [ns] $\\pm$ $k*{}$\\n{} == {}\".format(1/f_beacon/ns, beacon_time_delays/ns, None))\n",
    "\n",
    "for i, _ in enumerate(beacons):\n",
    "    l = ax.plot(time/ns, beacons[i], label=\"ch {}\".format(i), ls ='--', alpha=0.5)\n",
    "    \n",
    "    print(i, beacon_phase_delays[i], beacon_time_delays[i])\n",
    "\n",
    "    if not True:\n",
    "        corrected_beacon = sin_delay(f_beacon, time, phase=+beacon_phase_delays[i], t_delay=0)\n",
    "    else:\n",
    "        corrected_beacon = sin_delay(f_beacon, time, t_delay=beacon_time_delays[i], phase=0)\n",
    "        \n",
    "    ax.plot(time/ns, 2*beacon_amplitude*corrected_beacon, label='ch {} corrected'.format(i), color=l[0].get_color())\n",
    "    \n",
    "    # indicate start of uncorrected beacons\n",
    "    ax.axvline(beacon_time_delays[i]/ns, color=l[0].get_color())\n",
    "\n",
    "print(\"sum:\", np.sum(beacon_time_delays)/ns)\n",
    "print(\"diff:\", np.diff(beacon_time_delays)/ns)\n",
    "\n",
    "ax.legend(ncol=2)\n",
    "ax.margins(y=0.3)\n",
    "if True:\n",
    "    ax.set_xlim(time[0]/ns - 1, time[2*samplerate//f_beacon]/ns)\n",
    "\n",
    "fig.show()\n",
    "\n",
    "\n",
    "#    ax.plot((double_signal_time) * ns, signal_2(double_signal_time + calc_shift), 'r--', label='Recovered')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### 1.2 Impulse vs beacon delays\n",
    "\n",
    "Find the delay within a single beacon period"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "#def find_beacon_impulse_delay(samplerate, f_beacon, impulse, init_phase=0):\n",
    "def find_beacon_impulse_phase_delay(samplerate, f_beacon, reference_beacon, impulse, **lag_kwargs):\n",
    "    \"\"\"\n",
    "    Return phase delay of `beacon` with respect to `reference_beacon`.\n",
    "    Note that the returned value can be off by a multiple of $2\\pi$.\n",
    "    \n",
    "    Parameters\n",
    "    ==========\n",
    "    samplerate : float\n",
    "        Samplerate of both reference_beacon and delayed_beacon\n",
    "    f_beacon : float\n",
    "        Frequency of the beacons\n",
    "    reference_beacon : ndarray\n",
    "        The beacon to use as a reference\n",
    "    beacon : ndarray\n",
    "        The beacon to find the delay for\n",
    "    \"\"\"\n",
    "    \n",
    "    calc_lag, _ = find_best_lag(reference_beacon, impulse, **lag_kwargs)\n",
    "    \n",
    "    return 2*np.pi* f_beacon * calc_lag / samplerate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Beacon Impuls delays: $[ 867.07957239 1416.22996824]$ns\n"
     ]
    }
   ],
   "source": [
    "impulse_beacon_phase_delays = np.empty( len(impulses) )\n",
    "\n",
    "for i, _ in enumerate(impulses):\n",
    "    impulse_beacon_phase_delays[i] = find_beacon_impulse_phase_delay(\n",
    "        beacon_samplerate, f_beacon, \n",
    "        ref_beacon, impulses[i]\n",
    "    )\n",
    "\n",
    "print(\"Beacon Impuls delays: ${}$ns\".format(impulse_beacon_phase_delays/f_beacon/ns))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAVt0lEQVR4nO3df6zdd33f8ecLQxItUDDYk5DjEENdkhS1CRwZKqTBBEkc/rCRilqnQyQsnTVGYIOtUhCagsw/FDbRsaUlpvWASsOBVNpuO1CaEhhbh8HHIwuxK8PFUOKZKQaHaF2yZE7e++N8s3t8fa/v19fH9zj383xIRz7fz/fz+Z7P/eje8/L35ydVhSSpXc+bdgckSdNlEEhS4wwCSWqcQSBJjTMIJKlxBoEkNW7JIEiyJ8kjSR5aZH2SfCrJbJIHk7x2bN3NSb7fvW6eZMclSZPRZ4/gs8DWM6y/EdjcvXYCfwCQ5KXAHcDrgS3AHUnWnktnJUmTt2QQVNU3gBNnqLId+HyN7ANekuTlwA3AfVV1oqoeBe7jzIEiSZqC509gGxuAh8eWj3Zli5WfJslORnsTXHrppa+78sorJ9AtSWrHgQMHflpV65fTdhJBkAXK6gzlpxdW7QZ2AwwGgxoOhxPoliS1I8lfL7ftJK4aOgpsHFu+DDh2hnJJ0gVkEkEwA7yru3roDcBjVfUT4F7g+iRru5PE13dlkqQLyJKHhpJ8AXgzsC7JUUZXAr0AoKo+DXwZeBswCzwOvLtbdyLJR4H93aZ2VdWZTjpLkqZgySCoqpuWWF/AexdZtwfYs7yuSZJWgncWS1LjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIa1ysIkmxNcjjJbJLbF1j/ySQPdK/vJfn52Lqnx9bNTLLzkqRz12eqyjXAncB1jCak359kpqoOPVunqj4wVv99wLVjm3iiqq6ZXJclSZPUZ49gCzBbVUeq6ilgL7D9DPVvAr4wic5Jks6/PkGwAXh4bPloV3aaJK8ANgH3jxVfkmSYZF+Sty/SbmdXZ3j8+PGeXZckTUKfIMgCZbVI3R3APVX19FjZ5VU1AH4L+L0krzptY1W7q2pQVYP169f36JIkaVL6BMFRYOPY8mXAsUXq7mDeYaGqOtb9ewT4OqeeP5AkTVmfINgPbE6yKclFjL7sT7v6J8mrgbXAN8fK1ia5uHu/DngjcGh+W0nS9Cx51VBVnUxyG3AvsAbYU1UHk+wChlX1bCjcBOytqvHDRlcBdyV5hlHofGz8aiNJ0vTl1O/t6RsMBjUcDqfdDUl6TklyoDsfe9a8s1iSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1LheQZBka5LDSWaT3L7A+luSHE/yQPf67bF1Nyf5fve6eZKdlySduyWnqkyyBrgTuI7RRPb7k8wsMOXk3VV127y2LwXuAAZAAQe6to9OpPeSpHPWZ49gCzBbVUeq6ilgL7C95/ZvAO6rqhPdl/99wNbldVWSdD70CYINwMNjy0e7svl+PcmDSe5JsvFs2ibZmWSYZHj8+PGeXZckTUKfIMgCZfNnvP9T4Iqq+hXgL4DPnUVbqmp3VQ2qarB+/foeXZIkTUqfIDgKbBxbvgw4Nl6hqn5WVU92i58BXte3rSRpuvoEwX5gc5JNSS4CdgAz4xWSvHxscRvwV937e4Hrk6xNsha4viuTJF0glrxqqKpOJrmN0Rf4GmBPVR1MsgsYVtUM8P4k24CTwAnglq7tiSQfZRQmALuq6sR5+DkkScuUqtMO2U/VYDCo4XA47W5I0nNKkgNVNVhOW+8slqTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXG9giDJ1iSHk8wmuX2B9R9McqibvP6rSV4xtu7pJA90r5n5bSVJ07XkDGVJ1gB3AtcxmoN4f5KZqjo0Vu07wKCqHk/yHuDjwG92656oqmsm3G9J0oT02SPYAsxW1ZGqegrYC2wfr1BVX6uqx7vFfYwmqZckPQf0CYINwMNjy0e7ssXcCnxlbPmSJMMk+5K8faEGSXZ2dYbHjx/v0SVJ0qQseWgIyAJlC050nOSdwAB401jx5VV1LMkrgfuTfLeqfnDKxqp2A7thNGdxr55Lkiaizx7BUWDj2PJlwLH5lZK8FfgwsK2qnny2vKqOdf8eAb4OXHsO/ZUkTVifINgPbE6yKclFwA7glKt/klwL3MUoBB4ZK1+b5OLu/TrgjcD4SWZJ0pQteWioqk4muQ24F1gD7Kmqg0l2AcOqmgE+AbwQ+FISgB9X1TbgKuCuJM8wCp2PzbvaSJI0Zam6sA7JDwaDGg6H0+6GJD2nJDlQVYPltPXOYklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS43oFQZKtSQ4nmU1y+wLrL05yd7f+W0muGFv3oa78cJIbJtd1SdIkLBkESdYAdwI3AlcDNyW5el61W4FHq+oXgU8Cv9u1vZrRHMe/DGwFfr/bniTpAtFnj2ALMFtVR6rqKWAvsH1ene3A57r39wBvyWjy4u3A3qp6sqp+CMx225MkXSCWnLwe2AA8PLZ8FHj9YnW6ye4fA17Wle+b13bD/A9IshPY2S0+meShXr1f/dYBP512Jy4QjsUcx2KOYzHn1ctt2CcIskDZ/BnvF6vTpy1VtRvYDZBkuNwJmFcbx2KOYzHHsZjjWMxJMlxu2z6Hho4CG8eWLwOOLVYnyfOBFwMneraVJE1RnyDYD2xOsinJRYxO/s7MqzMD3Ny9fwdwf1VVV76ju6poE7AZ+PZkui5JmoQlDw11x/xvA+4F1gB7qupgkl3AsKpmgD8C/jjJLKM9gR1d24NJvggcAk4C762qp5f4yN3L/3FWHcdijmMxx7GY41jMWfZYZPQfd0lSq7yzWJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJalyfqSr3JHlkscliMvKpbl7iB5O8dmzdzUm+371uXqi9JGm6+uwRfJbRfMOLuZHR46U3M5pl7A8AkrwUuIPRbGZbgDuSrD2XzkqSJm/JIKiqbzB6tPRitgOfr5F9wEuSvBy4Abivqk5U1aPAfZw5UCRJU9BnqsqlLDSn8YYzlJ9mfM7iSy+99HVXXnnlBLolSe04cODAT6tq/XLaTiIIzmm+Yjh1zuLBYFDD4bKn3pSkJiX56+W2ncRVQ4vNS+x8xZL0HDCJIJgB3tVdPfQG4LGq+gmjqS2vT7K2O0l8fVcmSbqALHloKMkXgDcD65IcZXQl0AsAqurTwJeBtwGzwOPAu7t1J5J8FNjfbWpXVZ3ppLMkaQr6TF5/0xLrC3jvIuv2AHuW1zVJ0krwzmJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuN6BUGSrUkOJ5lNcvsC6z+Z5IHu9b0kPx9b9/TYuplJdl6SdO76TFW5BrgTuI7RhPT7k8xU1aFn61TVB8bqvw+4dmwTT1TVNZPrsiRpkvrsEWwBZqvqSFU9BewFtp+h/k3AFybROUnS+dcnCDYAD48tH+3KTpPkFcAm4P6x4kuSDJPsS/L2Rdrt7OoMjx8/3rPrkqRJ6BMEWaCsFqm7A7inqp4eK7u8qgbAbwG/l+RVp22sandVDapqsH79+h5dkiRNSp8gOApsHFu+DDi2SN0dzDssVFXHun+PAF/n1PMHkqQp6xME+4HNSTYluYjRl/1pV/8keTWwFvjmWNnaJBd379cBbwQOzW8rSZqeJa8aqqqTSW4D7gXWAHuq6mCSXcCwqp4NhZuAvVU1ftjoKuCuJM8wCp2PjV9tJEmavpz6vT19g8GghsPhtLshSc8pSQ5052PPmncWS1LjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1LheQZBka5LDSWaT3L7A+luSHE/yQPf67bF1Nyf5fve6eZKdlySduyVnKEuyBrgTuI7R/MX7k8wsMNPY3VV127y2LwXuAAaMJrw/0LV9dCK9lySdsz57BFuA2ao6UlVPAXuB7T23fwNwX1Wd6L787wO2Lq+rkqTzoU8QbAAeHls+2pXN9+tJHkxyT5KNZ9M2yc4kwyTD48eP9+y6JGkS+gRBFiibP9HxnwJXVNWvAH8BfO4s2lJVu6tqUFWD9evX9+iSJGlS+gTBUWDj2PJlwLHxClX1s6p6slv8DPC6vm0lSdPVJwj2A5uTbEpyEbADmBmvkOTlY4vbgL/q3t8LXJ9kbZK1wPVdmSTpArHkVUNVdTLJbYy+wNcAe6rqYJJdwLCqZoD3J9kGnAROALd0bU8k+SijMAHYVVUnzsPPIUlaplSddsh+qgaDQQ2Hw2l3Q5KeU5IcqKrBctp6Z7EkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXG9giDJ1iSHk8wmuX2B9R9McijJg0m+muQVY+ueTvJA95qZ31aSNF1LTlWZZA1wJ3Ado8no9yeZqapDY9W+Awyq6vEk7wE+Dvxmt+6Jqrpmwv2WJE1Inz2CLcBsVR2pqqeAvcD28QpV9bWqerxb3AdcNtluSpLOlz5BsAF4eGz5aFe2mFuBr4wtX5JkmGRfkrcv1CDJzq7O8Pjx4z26JEmalCUPDQFZoGzBGe+TvBMYAG8aK768qo4leSVwf5LvVtUPTtlY1W5gN4wmr+/Vc0nSRPTZIzgKbBxbvgw4Nr9SkrcCHwa2VdWTz5ZX1bHu3yPA14Frz6G/kqQJ6xME+4HNSTYluQjYAZxy9U+Sa4G7GIXAI2Pla5Nc3L1fB7wRGD/JLEmasiUPDVXVySS3AfcCa4A9VXUwyS5gWFUzwCeAFwJfSgLw46raBlwF3JXkGUah87F5VxtJkqYsVRfWIfnBYFDD4XDa3ZCk55QkB6pqsJy23lksSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWpcryBIsjXJ4SSzSW5fYP3FSe7u1n8ryRVj6z7UlR9OcsPkui5JmoQlgyDJGuBO4EbgauCmJFfPq3Yr8GhV/SLwSeB3u7ZXM5rj+JeBrcDvd9uTJF0g+uwRbAFmq+pIVT0F7AW2z6uzHfhc9/4e4C0ZTV68HdhbVU9W1Q+B2W57kqQLxJKT1wMbgIfHlo8Cr1+sTjfZ/WPAy7ryffPabpj/AUl2Aju7xSeTPNSr96vfOuCn0+7EBcKxmONYzHEs5rx6uQ37BEEWKJs/4/1idfq0pap2A7sBkgyXOwHzauNYzHEs5jgWcxyLOUmGy23b59DQUWDj2PJlwLHF6iR5PvBi4ETPtpKkKeoTBPuBzUk2JbmI0cnfmXl1ZoCbu/fvAO6vqurKd3RXFW0CNgPfnkzXJUmTsOShoe6Y/23AvcAaYE9VHUyyCxhW1QzwR8AfJ5lltCewo2t7MMkXgUPASeC9VfX0Eh+5e/k/zqrjWMxxLOY4FnMciznLHouM/uMuSWqVdxZLUuMMAklq3NSC4FweW7Ha9BiLDyY5lOTBJF9N8opp9HMlLDUWY/XekaSSrNpLB/uMRZLf6H43Dib5dyvdx5XS42/k8iRfS/Kd7u/kbdPo5/mWZE+SRxa71yojn+rG6cEkr+214apa8Rejk84/AF4JXAT8d+DqeXX+EfDp7v0O4O5p9PUCGYu/C/yt7v17Wh6Lrt6LgG8wullxMO1+T/H3YjPwHWBtt/y3p93vKY7FbuA93furgR9Nu9/naSz+DvBa4KFF1r8N+Aqje7jeAHyrz3antUdwLo+tWG2WHIuq+lpVPd4t7mN0P8Zq1Of3AuCjwMeB/7OSnVthfcbiHwB3VtWjAFX1yAr3caX0GYsCfqF7/2JW6f1KVfUNRldmLmY78Pka2Qe8JMnLl9rutIJgocdWzH/0xCmPrQCefWzFatNnLMbdyijxV6MlxyLJtcDGqvqzlezYFPT5vfgl4JeS/GWSfUm2rljvVlafsfgI8M4kR4EvA+9bma5dcM72+wTo94iJ8+FcHlux2vT+OZO8ExgAbzqvPZqeM45FkucxerrtLSvVoSnq83vxfEaHh97MaC/xPyd5TVX9/Dz3baX1GYubgM9W1b9M8muM7mt6TVU9c/67d0FZ1vfmtPYIzuWxFatNr8dwJHkr8GFgW1U9uUJ9W2lLjcWLgNcAX0/yI0bHQGdW6Qnjvn8j/6Gq/m+Nnu57mFEwrDZ9xuJW4IsAVfVN4BJGD6RrzbIe6zOtIDiXx1asNkuORXc45C5GIbBajwPDEmNRVY9V1bqquqKqrmB0vmRbVS37YVsXsD5/I/+e0YUEJFnH6FDRkRXt5croMxY/Bt4CkOQqRkFwfEV7eWGYAd7VXT30BuCxqvrJUo2mcmiozuGxFatNz7H4BPBC4Evd+fIfV9W2qXX6POk5Fk3oORb3AtcnOQQ8DfxOVf1ser0+P3qOxT8FPpPkA4wOhdyyGv/jmOQLjA4FruvOh9wBvACgqj7N6PzI2xjN/fI48O5e212FYyVJOgveWSxJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBoVUnysiQPdK//meR/jC3/1/PwebckOZ7kD5fZ/hNdP//ZpPsm9TWtZw1J50V3Q9U1AEk+AvxNVf2L8/yxd1fVbctpWFW/k+R/T7pD0tlwj0DNSPI33b9vTvKfknwxyfeSfCzJ30vy7STfTfKqrt76JH+SZH/3emOPz7glyb8ZW/6z7vPWJPlskoe6z/jA+ftJpbPjHoFa9avAVYweX3IE+MOq2pLkHzN6hPE/Af4V8Mmq+i9JLmf0iIOrlvl51wAbquo1AElecq4/gDQpBoFatf/Zh3El+QHw5135d+ke5Aa8Fbh6bD6kX0jyoqr6X8v4vCPAK5P8a+A/jn2eNHUGgVo1/ijvZ8aWn2Hu7+J5wK9V1RNnsd2TnHrI9RKAqno0ya8CNwDvBX4D+PvL6Lc0cZ4jkBb358D/Pwmc5JoebX4EXJPkeUk2Mppm8dnHRD+vqv4E+OeM5p2VLgjuEUiLez9wZ5IHGf2tfAP4h0u0+Uvgh4wOMT0E/LeufAPwb7tZ1gA+NPnuSsvjY6ilc5DkFmCw3MtHu218hJW5zFVakIeGpHPzBHDjudxQBrwT8F4CTY17BJLUOPcIJKlxBoEkNc4gkKTGGQSS1Lj/B5T7qtMWY+TaAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Make a figure showing the corrected beacons\n",
    "fig, axes = plt.subplots(len(impulses),1, sharex=True)\n",
    "axes[-1].set_xlabel(\"Time [us]\")\n",
    "ax.set_title(\"Beacon Impuls delays: ${}$ns\".format(impulse_beacon_phase_delays/f_beacon/ns))\n",
    "\n",
    "for i, beacon in enumerate(beacons):\n",
    "    ax.set_ylabel(\"Amplitude [au]\")\n",
    "    ax.plot(time, beacon, label=\"ch {}\".format(i))\n",
    "\n",
    "#ax.plot(beacon_time, corrected_beacon, label='phase corrected (overlaps ch 0)')\n",
    "\n",
    "ax.legend()\n",
    "\n",
    "fig.show()"
   ]
  },
  {
   "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
}