mirror of
https://gitlab.science.ru.nl/mthesis-edeboone/m-thesis-introduction.git
synced 2024-11-13 10:03:32 +01:00
563 lines
111 KiB
Text
563 lines
111 KiB
Text
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"### Discrete Sine + Delta correlation\n",
|
|
"\n",
|
|
"Bandpass a delta peak, and try to correlate it with an appropriate sine wave.\n",
|
|
"\\\n",
|
|
"Show that we can lift the degeneracy in the correlations (within a window)."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 1,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"import matplotlib.pyplot as plt\n",
|
|
"import numpy as np\n",
|
|
"import scipy.fft as ft\n",
|
|
"import scipy.signal as signal\n",
|
|
"\n",
|
|
"# monkey patch correlation_lags into signal if it does not exist\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": 2,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def plot_spectrum( ax, spectrum, freqs, plot_complex=False, plot_power=False, plot_amplitude=None):\n",
|
|
" \"\"\" Plot a signal's spectrum on an Axis object\"\"\"\n",
|
|
" plot_amplitude = plot_amplitude or (not plot_power and not plot_complex)\n",
|
|
" alpha = 1\n",
|
|
" \n",
|
|
" ax.set_title(\"Spectrum\")\n",
|
|
" ax.set_xlabel(\"f (Hz)\")\n",
|
|
" ylabel = \"\"\n",
|
|
" if plot_amplitude or plot_complex:\n",
|
|
" ylabel = \"Amplitude\"\n",
|
|
" if plot_power:\n",
|
|
" if ylabel:\n",
|
|
" ylabel += \"|\"\n",
|
|
" ylabel += \"Power\"\n",
|
|
" ax.set_ylabel(ylabel)\n",
|
|
"\n",
|
|
" if plot_complex:\n",
|
|
" alpha = 0.5\n",
|
|
" ax.plot(freqs, np.real(spectrum), '.-', label='Real', alpha=alpha)\n",
|
|
" ax.plot(freqs, np.imag(spectrum), '.-', label='Imag', alpha=alpha)\n",
|
|
"\n",
|
|
" if plot_power:\n",
|
|
" ax.plot(freqs, np.abs(spectrum)**2, '.-', label='Power', alpha=alpha)\n",
|
|
" \n",
|
|
" if plot_amplitude:\n",
|
|
" ax.plot(freqs, np.abs(spectrum), '.-', label='Abs', alpha=alpha)\n",
|
|
"\n",
|
|
" ax.legend()\n",
|
|
"\n",
|
|
" return ax\n",
|
|
"\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()*N_samples).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": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Generate a delta peak and bandpass it"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 3,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Size: 25; Offset: 0; PeakLength: 1\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"# Generate a delta peak and bandpass it\n",
|
|
"samplerate = 500 # MHz\n",
|
|
"band = (30, 80) # MHz\n",
|
|
"\n",
|
|
"us = 1e3 # ns\n",
|
|
"ns = 1/us # us\n",
|
|
"\n",
|
|
"## Delta peak properties\n",
|
|
"imp_timelength = 50 * ns# ns\n",
|
|
"imp_peaklength = 1\n",
|
|
"imp_time = np.arange(0, imp_timelength, 1/samplerate)\n",
|
|
"imp_offset = [1400]\n",
|
|
"\n",
|
|
"orig_imp, imp_offset = deltapeak(imp_timelength, samplerate, offset=imp_offset, peaklength=imp_peaklength)\n",
|
|
"\n",
|
|
"## Bandpass it\n",
|
|
"imp, (fft, fft_freqs) = fft_bandpass(orig_imp, band, samplerate)\n",
|
|
"imp /= np.max(imp)\n",
|
|
"\n",
|
|
"#\n",
|
|
"print(\"Size: {}; Offset: {}; PeakLength: {}\".format(orig_imp.size, imp_offset[0], len(imp_offset)))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 4,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stderr",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"<ipython-input-4-5b3a8514ce9c>:22: UserWarning: Matplotlib is currently using module://ipykernel.pylab.backend_inline, which is a non-GUI backend, so cannot show the figure.\n",
|
|
" fig.show()\n"
|
|
]
|
|
},
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA8gAAAEWCAYAAAC6+b50AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABtSklEQVR4nO3deXxU1d3H8c8vOwGysEPCEgTZIZCIqKhQ9x0VXFpt0aptbWttq90XWx9ba2u1trVuj1VbHxFQccNdEa1aJRD2fU/YIQkJSch2nj9mggETMiGTuTOT7/v1yiuZO3f5nTMzufc359xzzDmHiIiIiIiISHsX43UAIiIiIiIiIuFACbKIiIiIiIgISpBFREREREREACXIIiIiIiIiIoASZBERERERERFACbKIiIiIiIgIoARZpN0ws6+Y2ZshOM4kMyto6+OIiIiIiASbEmSRKGNmE83sIzMrMbN9ZvYfMzvBOfe0c+5sr+MTERGJFk2dc9vweJvM7My22r+IQJzXAYhI8JhZCvAK8C1gJpAAnAoc9DIuERGRaBOO51wzi3PO1Xh1fJFooBZkkehyPIBz7hnnXK1zrsI596ZzbomZTTezD+tXNLOzzWy1/1vvB83sfTO7wf/cdDP70Mz+ZGZFZrbRzM5rsO11ZrbSzErNbIOZfSP0RRUREfFUc+fc/5jZX/3n2VVmdkb9hmaWamb/a2bbzazQzP7HzGIbPH9jg/PsCjMbZ2b/AvoBL5tZmZn9yMwGmJkzs6+b2Rbg3cZudWrY8mxmd5jZLDP7t3//S83seDP7qZntMrOtZqYeZ9JuKUEWiS5rgFoze9LMzjOz9MZWMrNuwGzgp0BXYDVw8hGrnehf3g24B/hfMzP/c7uAC4EU4DrgPjMbF+zCiIiIhLHmzrknAhvwnUd/DTxvZl38zz0J1ACDgLHA2UD9l9TTgDuAr+I7z14M7HXOXQtsAS5yznVyzt3T4FinA8OAcwKM/SLgX0A6sAh4A19ekAH8Fng4wP2IRB0lyCJRxDm3H5gIOOBRYLeZvWRmPY9Y9XxguXPueX9XrAeAHUess9k596hzrhbfibw30NN/nFedc+udz/vAm/i6lYmIiLQLAZxzdwH3O+eqnXPP4vvS+QL/8+cBtzrnDjjndgH3AVf5t7sBuMc595n/PLvOObe5mXDu8O+rIsDwP3DOveG/BpgFdAfuds5VAzOAAWaWFuC+RKKKEmSRKOOcW+mcm+6cywRGAn2A+49YrQ+wtcE2Djhy5OkdDZ4v9//ZCcD/Tfkn/gFJivEl3N2CWQ4REZFw18w5t9B/fq232f98fyAe2G5mxf7z6MNAD/96fYH1LQxla/OrHGZng78rgD3+L8TrH4P/nC/S3ihBFolizrlVwBP4TtoNbQcy6x/4u05nEgAzSwSeA/4E9HTOpQFzATvadiIiItGskXNuRoNbk8B3//A2fMnsQaCbcy7N/5PinBvhX28rcFxThwlg+QEguf6B/97m7i0pi0h7pgRZJIqY2VAz+6GZZfof9wWuBj45YtVXgVFmNsXM4oBvA70CPEwCkAjsBmr8g3dpMA8REWlXAjjn9gBuMbN4/33Fw4C5zrnt+G5NutfMUswsxsyOM7PT/ds9BtxmZjnmM8jM+vuf2wkMbCa0NUCSmV1gZvHAL/Cdt0UkAEqQRaJLKb5BQf5rZgfwnaSXAT9suJJzbg8wDd/gW3uB4cACApiawjlXCtyCb0qLIuDLwEvBK4KIiEhEaO6c+19gMLAHuAuY6pzb63/uq/i+cF6B71w6G99YHzjnZvnX/z//MeYA9YN7/R74hb9r9m2NBeWcKwFuxpdoF+JrUT7yNioRaYIdfmuEiLRHZhaD7+T5Fefce17HIyIiEsnMbDpwg3NuotexiEjLqAVZpJ0ys3PMLM1/T/HP8N1DfGRXbBERERGRdkMJskj7dRK+UTL34JsPcUoLpocQEREREYk66mItIiIiIiIiglqQRURERERERACI8zqAUOrWrZsbMGCAt0GsXk1tbS2xw4d7G0cYOnDgAB07dvQ6jLCkummc6qVpqpvGBbte8vLy9jjnNL9oK6SlpblBgwZ5HUZQRMvnLlrKASpLuIqWskRLOSC6ytLac3O7SpAHDBjAggULvA1i0iSKi4tJ8zqOMDRv3jwmTZrkdRhhSXXTONVL01Q3jQt2vZjZ5qDtrJ3q2bOn9+fmIImWz120lANUlnAVLWWJlnJAdJWltedmdbEWERERERERoZ21IIeFX/yCzYsXk+Z1HCIiIiIiInIYJcihduaZFMWp2kVERERERMKNMrVQy8+n07p1ECV9/EUkNKqrqykoKKCysjKg9VNTU1m5cmUbRxV5jrVekpKSyMzMJD4+vg2iEhERkXChBDnUbr2VQcXFcMMNXkciIhGkoKCAzp07M2DAAMys2fVLS0vp3LlzCCKLLMdSL8459u7dS0FBAVlZWW0UmYiIiIQDDdIlIhIBKisr6dq1a0DJsQSXmdG1a9eAW+9FREQkcilBFhGJEEqOvaO6FxERaR+UIIuIiIiIiIigBFlERAIUGxtLdnY2Y8aMYdy4cXz00UdB2e+mTZsYOXJkUPYVDJMmTWLBggVehyEiIiIe0CBdofa737Fh4ULGeR2HiEgLdejQgfz8fADeeOMNfvrTn/L+++97G5SIiIhIEKkFOdROPpn9YdRSIiJyLPbv3096ejoAZWVlnHHGGYwbN45Ro0bx4osvAr6W4WHDhnHjjTcyYsQIzj77bCoqKgDIy8tjzJgxnHTSSfz9738/tN8nnniCSy65hHPPPZchQ4bwm9/85tBzU6ZMIScnhxEjRvDII48AUFtby/Tp0xk5ciSjRo3ivvvuA+CBBx5g+PDhjB49mquuugqAAwcOcPPNN3PCCScwduzYQ3FWVFRw1VVXMXr0aK688spDMYqIiEj7oxbkUPvoI1KWLdM8yCJyzH7z8nJWbNt/1HVqa2uJjY0NeJ/D+6Tw64tGHHWdiooKsrOzqaysZPv27bz77ruAb47gF154gZSUFPbs2cOECRO4+OKLAVi7di3PPPMMjz76KFdccQXPPfcc11xzDddddx1//etfOf3007n99tsPO86nn37KsmXLSE5O5oQTTuCCCy4gNzeXxx9/nC5dulBRUcEJJ5zA5ZdfzqZNmygsLGTZsmUAFBcXA3D33XezceNGEhMTDy276667OO200/jXv/5FcXEx48eP58wzz+Thhx8mOTmZJUuWsGTJEsaNUx8fERGR9kotyKH2s58x8LHHvI5CRKTF6rtYr1q1itdff52vfvWrOOdwzvGzn/2M0aNHc+aZZ1JYWMjOnTsByMrKIjs7G4CcnBw2bdpESUkJxcXFnH766QBce+21hx3nrLPOomvXrnTo0IHLLruMDz/8EPC1Co8ZM4YJEyawdetW1q5dy8CBA9mwYQPf/e53ef3110lJSQFg9OjRfOUrX+Hf//43cXG+74LffPNN7rvvPrKzs5k0aRKVlZVs2bKF+fPnc8011xzabvTo0W1elyIiIhKePG1BNrPHgQuBXc65L/Q7Nt+8Gn8BzgfKgenOuYX+5871PxcLPOacuztkgYuIeKi5ll6A0tJSOnfu3GYxnHTSSezZs4fdu3czd+5cdu/eTV5eHvHx8QwYMODQnMGJiYmHtomNjaWiogLn3FGnTTryOTNj3rx5vP3223z88cckJycfSnDT09NZvHgxb7zxBn//+9+ZOXMmjz/+OK+++irz58/npZde4s4772T58uU45/j3v//daAuxpnESERER8L4F+Qng3KM8fx4w2P9zE/APADOLBf7uf344cLWZDW/TSEVE5JBVq1ZRW1tL165dKSkpoUePHsTHx/Pee++xefPmo26blpZGamrqoZbhp59++rDn33rrLfbt20dFRQVz5szhlFNOoaSkhPT0dJKTk1m1ahWffPIJAHv27KGuro7LL7+cO++8k4ULF1JXV8fWrVuZPHky99xzD8XFxZSVlXHOOefw0EMP4ZwDYNGiRQCcdtpph2JYtmwZS5YsCWpdibSlvM1F/P29deRtLvI6lCaFe4zhHp+IhJanLcjOuflmNuAoq1wCPOV8VzOfmFmamfUGBgDrnHMbAMxshn/dFW0ccqvU1tWya38xMbV1XociItJi9fcgAzjnePLJJ4mNjeUrX/kKF110Ebm5uWRnZzN06NBm9/XPf/6T66+/nuTkZM4555zDnps4cSLXXnst69at48tf/jK5ubmMGjWKhx56iNGjRzNkyBAmTJgAQGFhIddddx11db7/q7///e+pra3lmmuuoaSkBOcc3//+90lLS+OXv/wl3/72txk9ejTOOQYMGMArr7zCt771La677jpGjx5NdnY248ePD27FibSRvM1FfPnRT6iurSMuJoafnDeEwT1b33Nk2Z5aYtfuDkKEsHZnKXe/tpqauuDGGKjmytIwvoS4GJ6+YQI5/dNDFp+IhB+r/ybdswB8CfIrTXSxfgW42zn3of/xO8CP8SXI5zrnbvAvvxY40Tn3nUb2cRO+1md69uyZM2PGjDYqSfPKqhzdbriFlATY9dhfPYsjXJWVldGpUyevwwhLqpvGtad6SU1NZdCgQQGv39JBusLF008/zcKFC7n33nvbZP+tqZd169ZRUlJy2LLJkyfnOedygxFbe9Lw3Ny9e/ecmTNnehxRcIT6f9Ir66uYvbY6ZMeLdjHAZYPjufC4BK9DaVQ0nfOipSzRUg6IrrK09twc7qNYN3ZTmDvK8i8udO4R4BGA3NxcN8nj0aPPv2A8MQfH8M+ck+jeObH5DdqRefPm4fXrE65UN41rT/WycuXKFt1T3Nb3ILeVpKQkEhIS2iz21tRLUlISY8eODXJE7VPDc/OQIUM8PzcHS6j/J1V338nstQswID42hl9fNJwhvVr/2Vm0aFHQ3uurd5Tym5dXHGpBDlaMgWquLKt3lHLHy8uprnXExBhXn3lC2LYgR9M5L1rKEi3lgOgqS2uFe4JcAPRt8DgT2AYkNLE87JXnlLJp3UCe/GgTt50zxOtwRETCyvTp05k+fbrXYYhEhDW7SgG44dQszh3ZO2iJXdmmWHIHdAnKvnIHdGFo7xQ+2bCXCQO7hjz5bK4s9fF95/8WEh8Tw7h+aaELTkTCkteDdDXnJeCr5jMBKHHObQc+AwabWZaZJQBX+dcNe5dsP8AVuz/hX59s5sDBGq/DERERkQjknOO5vALGD+jCzy8YHratngA5/dP59uRBYRtjfXxbispZVnj0OeZFJPp5miCb2TPAx8AQMysws6+b2TfN7Jv+VeYCG4B1wKPAzQDOuRrgO8AbwEpgpnNuecgLcAwuf34l3/jwWUoqqnn2s61ehyMiIiIRaOGWIjbsOcDU3EyvQ4kKF43pQ2JcDLPydG0m0t55PYr11c0874BvN/HcXHwJdERJiuuAxVYyvE8C//vhRr56Un/iYsO9IV9ERETCyey8ApITYrlgVG+vQ4kKqR3iOWdEL17M38bPLxhGYlzkDXIoIsGhzCzEOsQnAXDy0HIKiyt4del2jyMSERGRSFJRVcvLi7dz3sjedEwM9+FkIsfUnExKKqp5e8Uur0MREQ8pQQ6xpLgOACR32sBx3TvyyPwNeD3VlohIIAoKCrjkkksYPHgwxx13HN/73veoqqr6wnrbtm1j6tSpze7v/PPPp7i4+JhiueOOO/jTn/50TNuKRLrXl2+n7GAN09S9OqhOGdSN3qlJ6mYt0s4pQQ6xuJg4DGNj8QZuPHUgy7ft56P1e70OS0TkqJxzXHbZZUyZMoW1a9eyZs0aysrK+PnPf37YejU1NfTp04fZs2c3u8+5c+eSlpbWRhGLRK/ZeQX065LM+CCNNC0+sTHG5eMymb9mNzv3V3odjoh4RAlyiNnDD/PBN65hY/FGpozNoFunRB56f73XYYmIHNW7775LUlIS1113HQCxsbHcd999PP744zz44INMmzaNiy66iLPPPptNmzYxcuRIAMrLy7niiisYPXo0V155JSeeeCILFiwAYMCAAezZs4dNmzYxbNgwbrzxRkaMGMHZZ59NRUUFAI8++ignnHACY8aM4fLLL6e8vNybChAJEwVF5Xy0fi+Xj8skJsa8DifqXJ6TSZ2D5xcWeh2KiHhEN66E2pAhlK/oy4ZtM0mKj+W6UwbwxzdWs2Lbfob3SfE6OhGJALe+fiv5O/KPuk5tbS2xsYEPMpPdK5v7z72/yeeXL19OTk7OYctSUlLo168fNTU1fPzxxyxZsoQuXbqwadOmQ+s8+OCDpKens2TJEpYtW0Z2dnaj+1+7di3PPPMMjz76KFdccQXPPfcc11xzDZdddhk33ngjAL/4xS/43//9X7773e8GXC6RaPP8wkKcg8tzMrwOJSpldevICQPSmZ23lW+ePhAzfQkh0t6oBTnUXn6ZgXnL2Vi8Eecc15zYn+SEWB79YIPXkYmINMk51+iFYv3ys846iy5dvtjd88MPP+Sqq64CYOTIkYwePbrR/WdlZR1KnnNycg4l2cuWLePUU09l1KhRPP300yxfHhEz+om0ibo6x+y8Ak4+riuZ6clehxO1puZksn73ARZtLfY6FBHxgFqQQ+3eexmxbzeVl1eyo2wHvTv35qoT+vHUx5u4/Zwh9Enr4HWEIhLmjtbSW6+0tJTOnTsH7ZgjRozgueeeO2zZ/v372bp1K7GxsXTs2LHR7QIdhDAxMfHQ37GxsYe6WE+fPp05c+YwZswYnnjiCebNm3dsBRCJAp9u2seWfeV8/6zBXocS1S4Y3Yc7XlrBrAUFjOuX7nU4IhJiakH2QIy/2jcU+VqNr584AAc8/uFGD6MSEWnaGWecQXl5OU899RTg68L9wx/+kOnTp5Oc3HRL1sSJE5k5cyYAK1asYOnSpS06bmlpKb1796a6upqnn3762AsgEgVm5xXQKTGOc0do7uO21CkxjvNG9eKVxduorK71OhwRCTElyB6I8XdT3FjsS4gz05O5cHRvnvl0CyUV1V6GJiLSKDPjhRdeYNasWQwePJjjjz+epKQkfve73x11u5tvvpndu3czevRo/vCHPzB69GhSU1MDPu6dd97JiSeeyFlnncXQoUNbWwyRiHXgYA1zl27nwtG96ZAQ+PgCcmym5mRSerCGN5bv8DoUEQkxdbH2gNnhLcgAN502kBfzt/H0fzdz86RBXoUmItKkvn378vLLL39h+fTp05k+ffqhxwMGDGDZsmUAJCUl8e9//5ukpCTWr1/PGWecQf/+/QEO3WfcrVu3Q+sD3HbbbYf+/ta3vsW3vvWtLxzzjjvuCEKJRCLHq0u3U15Vq7mPQ2RCVlcy0zswa0EBl2RrQDSR9kQJsgcMyOiccagFGWBEn1ROHdyNf/5nE1+fmEVinL4dFpHIV15ezuTJk6mursY5xz/+8Q8SEhK8Dksk4szOK2Bgt466JzZEYvxzIj/w7loKiyvI0BgxIu2GuliH2r/+xcqf/Yys9KzDWpDB14q8u/QgLy7a5lFwIiLB1blzZxYsWMDixYtZsmQJ5513ntchiUSczXsP8OnGfVyek6lph0Joak4mzsHzeQVehyIiIaQEOdT69uVgjx4MTB/IxqLDB+WaOKgbw3un8MgHG6irC2zkVxFpPwIdEVqCT3UvXnour4AYg8vHqXt1KPXtksxJA7sye2GB/geItCNKkEPt2Wfp/u67ZKVlUbC/gIM1Bw89ZWbcdNpA1u0q491VuzwMUkTCTVJSEnv37tVFmgecc+zdu5ekpCSvQ5F2qK7O8dzCQiYO7k6vVL0HQ21qTiab95bz2aYir0MRkRDRPcih9o9/kFFczMDLf4DDsaVkC4O7fj6f4QWje3PP66t4ZP4Gzhze08NARSScZGZmUlBQwO7duwNav7KyUgldI461XpKSksjMVOudhN5H6/dSWFzBT87TKO5eOG9UL3790nJmLdjK+KwuXocjIiGgBNkjWWlZgG8k64YJcnxsDNdPzOJ/Xl3Joi1FjNVgHCICxMfHk5WVFfD68+bNY+zYsW0YUWRSvUikmZ23lZSkOM7Sl+aeSE6I44JRvXl5yTbuuHgEHRN16SwS7dTF2iMD0wcCfGGgLoCrxvcjJSmOR+Z/8TkRERFpH/ZXVvPash1cnN2HpHjNbuGVqbmZlFfV8toyzYks0h4oQfZI7869SYxNPGyqp3qdEuO4ZkJ/Xl++g017DngQnYiIiHjt1SXbOVhTx7Scvl6H0q7l9k9nQNdkZudt9ToUEQkBTxNkMzvXzFab2Toz+0kjz99uZvn+n2VmVmtmXfzPbTKzpf7nFoQ++taJsRgGpA1otAUZYPrJA4iPieHRD9SKLCIi0h7NWrCVwT06MToz1etQ2jUzY2pOJp9s2MeWveVehyMibcyzBNnMYoG/A+cBw4GrzWx4w3Wcc390zmU757KBnwLvO+f2NVhlsv/53FDF3WqzZ7P8N78BfN2sG2tBBuiRksSlYzOYnVfAnrKDja4jIiIi0WndrjIWbilmWq7mPg4Hl43LxAxmL9ScyCLRzssW5PHAOufcBudcFTADuOQo618NPBOSyNpSt25Up/q+Cc5Ky2qyBRngxtOyOFhTx1Mfbw5VdCIiIhIGnltYQGyMMWVshtehCNAnrQMTB3XjubwC6uo03Z5INDOv5tQ0s6nAuc65G/yPrwVOdM59p5F1k4ECYFB9C7KZbQSKAAc87Jx7pInj3ATcBNCzZ8+cGTNmtEVxAtbr9deprKykeMoUdh7YScH+ArJ7ZRNrjQ++cX9eJeuKa7l3UjKJsdH9DXJZWRmdOnXyOoywpLppnOqlaaqbxgW7XiZPnpwXUb2YwkTDc3P37t1zZs6c6XFEwRGs91edc/xgXgX9U2L4fk7op2uLpv8fwSzLJ9tqeGjJQX50QhLDu4Z+0DS9LuEnWsoB0VWW1p6bvRyrvrFsr6ls/SLgP0d0rz7FObfNzHoAb5nZKufc/C/s0Jc4PwKQm5vrJk2a1MqwW+mOOyguLibt/vt5fuXz3DbzNvIm5TGu97hGV+84YB/THvqYnclZfPWkAaGNNcTmzZuH569PmFLdNE710jTVTeNUL+Gh4bl5yJAh3p+bgyRY7695q3dRfPAz7j5nDJNG9m59YC09fhR9ToJZlgnVtTy95m3W1XTl5knZQdlnS+h1CT/RUg6IrrK0lpddrAuAhsMyZgLbmlj3Ko7oXu2c2+b/vQt4AV+X7YhSP9XTxqLG70MG38iJY/ul8egHG6iprQtVaCIiIuKRWXkFpCfH86Whmvs4nCTFx3LRmD7MXbad0spqr8MRkTbiZYL8GTDYzLLMLAFfEvzSkSuZWSpwOvBig2Udzaxz/d/A2cCykEQdRFlpWUDjcyHXMzO+cdpAtu6r4PXlmn9PROSZT7fw6pLtXoch0iaKy6t4a/lOLsnOICFOs3GGm2k5mVRW1+l/kEgU8+w/r3OuBvgO8AawEpjpnFtuZt80s282WPVS4E3nXMMJgXsCH5rZYuBT4FXn3Ouhij1YUpNS6dKhS5MjWdc7a3gveqcm6Z+xiAjwj3nreXVpUx2ORCLby4u3UVVbx7TcTK9DkUZk901jUI9OzM7TaNYi0crLe5Bxzs0F5h6x7KEjHj8BPHHEsg3AmDYOLySaG8kaIDbGGNcvncUFxaEJSkQkTJWUV7NlXzlXje/b/MoiEWhWXgHDeqcwoo/mPg5H9XMi3/3aKjbsLmNg9+gY1EhEPqe+O6E2dy5L7r770MOjzYXc0MiMVAqKKig6UNWW0YmIhLVl20oAGJWh5EGiz+odpSwpKGFajlqPw9llYzOIjTG1IotEKSXIoZacTF3S51M2ZKVlsal4E7V1tUfdrP5isP7iUESkPVpaqARZotfsvK3Ex2ru43DXIyWJ04/vzvMLC6nVnMgiUUcJcqg9+CB95sw59HBg+kCqaqvYVnr0++nqLwbrLw5FRNqjpQUl9O3SgbTkBK9DEQmq6to6XlhUyJeG9qBLR72/w93UnEx27K/kw3V7vA5FRIJMCXKozZxJj3nzDj08NNVTM92sU5Pj6dclmWVKkEWkHVtaWKLWY4lK76/ezZ6yKqbl6P76SHDGsB6kJcerm7VIFFKC7LGs9Oaneqo3KiNVLcgi0m7VD9A1UgmyRKFZeVvp1imB04d09zoUCUBiXCyXjOnDG8t3UFKuOZFFookSZI/1S+1HjMWwsSiwgbq27quguFwDdYlI+6MBuiRa7S07yDsrd3Hp2AziY3VpFimm5falqqaOl5Zo2jmRaKL/wh5LiE0gMyWTDcWBtSCD7kMWkfap/n/fSE1/I1Hmxfxt1NQ5pqp7dUQZ0SeFob06q5u1SJRRghwGBqYPDLAFOQVQgiwi7dPSwhIy0zuQrgGMJMrMyitgdGYqQ3p19joUaYH6OZEXby1m7c5Sr8MRkSBRghxq8+aRf//9hy3KSssK6B7ktOQE+nbpoIG6RKRdWqYBuiQKLd9Wwsrt+zX3cYSaMjaDOM2JLBJVlCCHgYHpA9letp2K6opm19VAXSLSHpWUV7N5bzmjMpUgS3SZtaCAhNgYLhrTx+tQ5Bh065TI5KE9eH5RITW1dV6HIyJBoAQ51P70J/o+++xhi7LSfCNZbyre1OzmozLSNFCXiLQ7GqBLolFVTR0v5hdy1oiemts7gk3LyWR36UHeX7Pb61BEJAiUIIfaK6/Q9eOPD1sU6FzI8PnF4bLC/cGPTUQkTGmALolG767aSVF5tbpXR7jJQ3vQtWOCulmLRAklyGGgPkEO5D5kDdQlIu2RBuiSaDRrQQE9UxI5dbDmPo5k8bExTBmbwdsrd7LvgHr4iUQ6JchhoEfHHiTHJ2ugLhGRJmiALok2u0ormbdmN5eNyyQ2xrwOR1ppWm4m1bWOF/MLvQ5FRFpJCXIYMDOy0rIC6mINvm7WSwqL2zYoEZEwUVLhG6BrpBJkiSJzFhVSW+eYqu7VUWForxRGZaSqm7VIFFCCHGodOlCbmPiFxVnpgU31BDAyI1UDdYlIu7G8UAN0SXRxzjFrQQHj+qVxXPdOXocjQTI1J5Pl2/azYpvGiRGJZEqQQ+2111j6hz98YfHAtIFsLNqIc67ZXWigLhFpT5YoQZYos6SghLW7ypiW29frUCSILh7Th4TYGLUii0Q4JchhIis9i9KqUvZW7G123fpRXDVQl4i0BxqgS6LNrLytJMXHcMHo3l6HIkGU3jGBM4f3YE5+IVU1mhNZJFJ5miCb2blmttrM1pnZTxp5fpKZlZhZvv/nV4FuG7buvJP+Tz31hcWHpnoqav4+5PSOGqhLRNoPDdAl0aSyupaX8rdx7ohepCTFex2OBNm0nL7sO1DFu6t2eR2KiBwjzxJkM4sF/g6cBwwHrjaz4Y2s+oFzLtv/89sWbht+3nmH9IULv7A4Ky0LCGyqJ/B1NVQLsohEOw3QJdHmrRU72V9Zo+7VUerUwd3o0TlR3axFIpiXLcjjgXXOuQ3OuSpgBnBJCLYNS1npvgQ50JGsR2aksmVfOSXl1W0ZloiIpzRAl0SbWXkFZKR14KSBXb0ORdpAXGwMl47L4L3Vu9hdetDrcETkGMR5eOwMYGuDxwXAiY2sd5KZLQa2Abc555a3YFvM7CbgJoCePXsyb9681kfeCtnFxdTW1jYax/1D7ydtb1pAMdbtqQXg33PnM6JbbJCj9EZZWZnnr0+4Ut00TvXStGipm7kbfaP179+0jHnbWj9XbLTUS6RreG7u3r171Lwmzb2/iirr+GBNBRcdF8/8+e+HLrAWiqbPiRdl6V9bR22d497Z8zk3K3jd6PW6hJ9oKQdEV1lay8sEubErnSOHcF4I9HfOlZnZ+cAcYHCA2/oWOvcI8AhAbm6umzRp0rHGGxxpaRQXF9NYHD957Cd0rurMW5e81exuxhyo4k8L3iKm2wAmTTquDQINvXnz5jVaL6K6aYrqpWnRUjezty0kI62YC8+eHJT9RUu9RLqG5+YhQ4Z4f24OkubeX39/bx2O1fzwslPo37Vj6AJroWj6nHhVlpmb/8Oi4lp+f/qpmLX+yz3Q6xKOoqUcEF1laS0vu1gXAA1vwMnE10p8iHNuv3OuzP/3XCDezLoFsm3Y6tqV6pSURp8amD4w4HuQ0zsmkJmugbpEJLot1QBdEiWcc8zOK2B8VpewTo4lOKblZrJ6Z6nGixGJQF4myJ8Bg80sy8wSgKuAlxquYGa9zP+1m5mNxxfv3kC2DVvPPcfy3/620aey0rLYXLyZmrqagHalgbpEJJrVD9A1KlMJskS+hVuK2LjnANNyMr0ORULgwtF9SIzTnMgikcizBNk5VwN8B3gDWAnMdM4tN7Nvmtk3/atNBZb570F+ALjK+TS6behLEVwD0wdS62op2B/YP1MN1CUi0UwDdEk0mbWggOSEWM4fpbmP24PUDvGcM6IXL+Zvo7K61utwRKQFvLwHub7b9Nwjlj3U4O+/AX8LdNuI8NOfkrVlCzTSx79+JOsNRRsYkDag2V2N9reqLNtWwimDugUzShERzy1VgixRoqKqlleWbOf8Ub3pmOjppZeE0LTcTF5avI13Vu7igtH6YkQkUnjZxbp9+vhjUpc33tg9MH0gABuLApzqqY/volHdrEUkGi0tLCEjrQPpHRO8DkWkVV5fvp2ygzXqXt3OnHxcN/qkJjErb2vzK4tI2FCCHEYyUzKJtdgWD9SlBFlEotEyDdAlUWLWggL6dUlmfFYXr0OREIqNMS4bl8n8NbvZUVLpdTgiEiAlyGEkLiaO/mn92VgcWAsy+AfqKlCCLCLRZX9lNZs0QJdEgYKicj5av5epOZlBm+5HIsfUnEzqHLywqNDrUEQkQEqQw0xLpnoCDdQlItGpfgq7kWpBlgj3XF4hZnDZuAyvQxEPDOjWkRMGpDMrbyvOOa/DEZEAKEEOtcxMDnbv3uTTWWlZLW5BBt9AXSIi0aK+Z4y6WEskq6tzzF64lZOP60pmerLX4YhHpuX0ZcPuAyzcUux1KCISACXIofbvf7Py5z9v8umB6QPZdWAXZVVlAe2u/uJR9yGLSDSpH6Criwbokgj26aZ9bN1XwbScvl6HIh46f3RvOsTHak5kkQihBDnMZKX5pnraVLwpoPXTOyaQkaaBukQkumiALokGsxYU0DkxjnNG9PI6FPFQp8Q4zhvVi1cWb6OiSnMii4Q7JcihduutDPpbo1M7A59P9dSS+5BHZ6Yeul9PRCTSaYAuiQYHDtbw2rLtXDimNx0SYr0ORzw2LacvpQdreHPFDq9DEZFmKEEOtfx8Oq1b1+TTWem+FuSWDtS1ea8G6hKR6KABuiQavLp0O+VVtUzV3McCnJjVhcz0DsxaoG7WIuGu2QTZzJLN7Jdm9qj/8WAzu7DtQ2ufunboSueEzmws0kBdItI+1SfI6mIdeXTN8LnZCwoY2L0j4/qlex2KhIGYGGNqTib/Wb+HwuIKr8MRkaMIpAX5n8BB4CT/4wLgf9osonbOzMhKz2JDceAtyBqoS0SiydLC/RqgK3LpmgHYtOcAn27ap7mP5TCXj8vEOXheg3WJhLVAEuTjnHP3ANUAzrkKQP/t29DA9IEtakHWQF0iEk2WFhQzMiPF6zDk2OiaAXhuYQExBpeNVfdq+VzfLsmcNLArsxcWaE5kkTAWSIJcZWYdAAdgZsfh+3ZYjsXxx1OeefQTZv1cyC355zkqQwN1iUjkOzRAl7pXR6p2f81QV+d4Lq+AUwd3p1dqktfhSJiZlpvJ5r3lfLapyOtQRKQJgSTIvwZeB/qa2dPAO8CP2jSqaPbII6y57bajrjIwfSDl1eXsOrAr4N2OyvQP1FWhgbpEJHJpgK6I1+6vGT5av5dtJZVMy1XrsXzRuSN70SkxjlkLtnodiog0Ia65FZxzb5nZQmACvm5S33PO7WnzyNqxhlM99ezUM6Bt6i8mlxeWcPKgbm0Wm0S/uUu3s3hr8TFte/aInuT07xLcgKRd0QBdkU3XDDArbyspSXGcOSyw87e0L8kJcVwwqjcvL9nGHRePoGNis5fiIhJiTX4qzWzcEYu2+3/3M7N+zrmFbRdWFLvpJo7ftg0mTWpylaw031RPG4s3clLfk5pcr6GGA3UpQZZjNW/1Lm5+eiEJsTHEtHASuJpax9P/3cLcW06lX9fktglQol79AF1dOyV6HYq0gK4ZfMqrHa8v28EVuX1Jitfcx9K4abmZPLtgK68t26FpwETC0NG+trrX/zsJyAUW4/s2eDTwX2Bi24YWpdasIbm4+KirDEgbALRsLuQuGqhLWmlP2UFum7WEIT078+J3TmnxxV1hcQXn3T+fW2YsYtY3TyI+VtOsS8stKyzRAF2RSdcMwKc7ajhYU6ekR44qp386Wd06MmvBVr1XRMJQk1ewzrnJzrnJwGZgnHMu1zmXA4wF1oUqwPaoQ3wHenfq3aKRrMHXiqwEWY6Fc47bZy1mf2U1f7k6+5haPjLSOvC7y0aRv7WYB95Z2wZRSrTbX1nNxj0H1L06AumawefDwhqO79mJ0Zl6D0vTzHxzIv934z627C33OhwROUIgTTxDnXNL6x8455YB2cE4uJmda2arzWydmf2kkee/YmZL/D8fmdmYBs9tMrOlZpZvZguCEU84GZg+sEVzIYMG6pJj99THm3lv9W5+dt5QhvY69ta7C0f3YWpOJn9/bx2fbtwXxAilPdAAXVGhza4Zwt26XWWsK65jWk5fzX0szbpsXAZmMHuh5kQWCTeBJMgrzewxM5tkZqeb2aPAytYe2Mxigb8D5wHDgavNbPgRq20ETnfOjQbuBB454vnJzrls51xua+MJN1npWS1uQW44UJdIoFbvKOWuuSuZPKQ7Xzt5QKv3d8fFI+jXJZlbZyyipFxf1kjgNEBXVGiTa4ZIMDvPN/fxJWP7eB2KRIDeqR2YOKgbz+UVUFenOZFFwkkgCfJ1wHLge8CtwAr/stYaD6xzzm1wzlUBM4BLGq7gnPvIOVc/UdwnQOTfqJGdTdmgQc2uNjBtIFv3b6WqtirgXTccqEskEJXVtdzyzCJSkuL447QxQWn16JQYx1+uGsuu0oP8bM7SFs3nLe3b0sL99ElN0gBdka2trhnCWm2d44VFBYzuFkuPzpr7WAIzLbcvhcUVfLJhr9ehiEgD5tXFq5lNBc51zt3gf3wtcKJz7jtNrH8bvq5b9etvBIoABzzsnDuydbl+u5uAmwB69uyZM2PGjKCXpaXKysro1KnTUdfZW7GXTcWbGNljJImxgV8s/nBeOcelxXBzduSdoAOpl/aqrerm3ysO8vaWGn6Qk8jo7sGdauKV9VXMXlvN10cmcGpmfFD3XU/vmaZFYt38eH45GZ1iuGVc2/3/Cna9TJ48OS8aezG1tYbn5u7du+fMnDnT44haZ8nuGv6cd5AbhjomDoisz11jIvH/R1PCuSxVtY7vvVdOdo9YvjG6+f974VyWloqWskRLOSC6ytLac3OzV8T+RPQLWbRzbuCxHrR+140sazRbN7PJwNc5fBTMU5xz28ysB/CWma1yzs1vJM5H8HfNzs3NdZOOMr1SqMybN4/m4pi/eT5Tn5jKm+PfZNJxR1+3odytC1i9o7TZ/YejQOqlvWqLunlv1S7efv0zrjtlALdcNCKo+wY49TTH1kc/4Zk1JVx77skM6NYx6MfQe6ZpkVY3+yur2fn6m1w78TgmTRrcZseJtHqJNIFeMzQ8Nw8ZMiQszs2tMevphXTpuJcJ/eKi4v0VTZ+TcC/LZfuX8tzCAsZNOIWUpKN/mRzuZWmJaClLtJQDoqssrRVIF+tc4AT/z6nAA8C/g3DsAqBvg8eZwLYjVzKz0cBjwCXOuUN9UJxz2/y/dwEv4OuyHf6uuYZhd93V7GoD033XEi2Z6glgdGYamzRQlzRjd+lBbp+9mKG9OvPjc4e2yTFiY4z7rswmPjaG781YRHVtXZscR6LD8sL9gAboigJtdc0QtorLq3hrxU4uye5DXIwG55KWmZbbl8rqOuYu2d78yiISEs0myM65vQ1+Cp1z9wNfCsKxPwMGm1mWmSUAVwEvNVzBzPoBzwPXOufWNFje0cw61/8NnA0sC0JMba+ggMTdu5tdrU/nPiTEJrCxWAN1SXDV1Tlum7WY0soaHrh67DFN6RSoPmkduPuyUSwuKOG+t9Y0v4G0WxqgKzq04TVD2Hpp8TaqajX3sRybMZmpDOrRiVl5Gs1aJFwE0sV6XIOHMfi+He7c2gM752rM7DvAG0As8LhzbrmZfdP//EPAr4CuwIP+wYNq/P3JewIv+JfFAf/nnHu9tTGFkxiLYUDagBa3IDccqOvkQd3aIjSJcE98tIn31+zmzktGcHzPVn+Um3XeqN5cmduXf7y/nlMHd+ek47q2+TEl8iwpLNEAXVGgra4ZwtnsvAKG905hRJ9U5ul7QGkhM2NaTia/f20VG3aXMbB7dNwDKhLJAhmV594Gf9fgm3rpimAc3Dk3F5h7xLKHGvx9A3BDI9ttAMYcuTzaZKVltbgFuUvHBDLSOmgka2nUyu37ufu1VZw5rAfXTOgfsuP+6qLhfLppHz+Ymc9r3zuVtOSEkB1bIsOywhJ1r44ObXbNEI5W7yhlSUEJv77oyFkqRQJ36dgM7nljNbPzCvhRG932JCKBC+Qe5K875yb7f85yzt0EBD73kByzgekDW9yCDDAyI+VQd0WRevVTOqUmx/OHy0cHZUqnQHVMjOOBq8ayp+wgP31eUz/J4fZXVrNxzwF1r44O7eqaYdaCrcTHGpdkZ3gdikSwHilJnH58d55fWEit5kQW8VwgCfLsAJdJIE46iZIRgY0YnJWWxb6KfZRUtizZHZWRyqa95eyv1EBd8rm7Xl3J2l1l3DttjCfdWEdlpvLDs4fw2rIdzFywNeTHl/B1aICuTCXIUaDdXDNU19YxJ7+QM4b2pEtH9YqR1pmWk8mO/ZV8uG6P16GItHtNdrE2s6HACCDVzC5r8FQKEHmT7IaL3/+ejfPmEUjn1vqRrDcWbyS7V3bAh6jvprissISTj9N9yAJvr9jJvz7ZzA0Tszjt+O6exXHTqQOZv2Y3d7y0gtwBXThO91oJGqArGrTHa4Z5q3ezp6xKg3NJUHxpWA/SkuOZtWArp3t4nhaRo7cgDwEuBNKAixr8jANubPPIhKz0LAA2FrXsPuRRDRJkkV37K/nRc0sY3juF288d4mksMTHGn6/IJjE+hltn5FNVo6mfxDeoYJ/UJLppgK5I1u6uGWbnbaVbp0QmDVEyI62XGBfLlOwM3lyxk5Jy9QAU8VKTLcjOuReBF83sJOfcxyGMKbpdfjkjdu+G+fObXfVY50Lu2inRP1DX/mMKUaJHXZ3jh7MWU15VwwNXZ5MY13ZTOgWqV2oSd182mm/+O49731rNT88b5nVI4jEN0BX52ts1w96yg7yzchfXT8wiLjaQu9VEmjc1J5MnPtrES0u2cW0IB9IUkcMdrYv1j5xz9wBfNrOrj3zeOXdLm0YWrfbuJX5/YIlrWlIa6UnpxzxQ19KC4hZvJ9Hl8f9s5IO1e7jr0pEM6hE+M62cO7IXV4/vxyPzN3D64O6akqwdK62sZsOeA1w6VoMcRbL2ds0wJ38bNXVO3aslqEb0SWFor87MXrBVCbKIh472tedK/+8FQF4jPxICWektn+oJNFCX+Frl/vD6Ks4a3pMvj+/ndThf8MsLh5HVrSPfn5lP0YGoHeRWmrFMA3RFi3Z1zTA7r4AxmakhmUte2g8zY1puXxYXlLBmZ6nX4Yi0W00myM65l/2/n2zsJ3Qhtm/HPtWT7kNuzyqqavnejEWkJyeEfEqnQCUn+KZ+2negih8/t0RTP7VTGqArOrSna4ZlhSWs3L5frcfSJqZk9yEuxpidV+B1KCLt1tG6WL8MNHnF6py7uE0iksNkpWXx8uqXqXN1xFjg9zmN0kjWR7Vqx34eeGctl43N5IxhPcIygWyN381dyfrdB/j3108M6+lHRmakcvs5Q/jd3FXMXLCVK08Iv5bu1qiqqeO5hQV8tH4vv714BOlh/Fp4ZWlhCb01QFfEa0/XDLPzCkiIi+HiMbotQIKva6dEvjS0B88vLORH5wzRPe4iHmgyQQb+FLIo2pMzzqBo40bSAlx9YPpADtYeZHvpdjJSAj8Zd+2USJ/UJA3U1YjyqhpufnohG3YfYO7SHQzrncJ3vzSIc0f0IiYm8hPld1d9PqXTxMHh/+XIDRMHMm/1bn7z8grGZ3Ulq1tHr0NqtcrqWp79bCsPvb+e7SWVh5Y9cm1O1H0Z01oaoCtqtItrhqqaOl7ML+Ts4T1JTY73OhyJUtNy+/Lmip28v2Y3Zwzr6XU4Iu3O0bpYv1//A3wMFAH7gI/9y+RY/PKXbP7qVwNePSvNP9XTMdyHPDIjVV2sG3HnKyvYuOcAT14/nnunjeFgdS03P72Qs++fz5xFhdTURu7UQ3vKDvKj2UsY2quz51M6BSomxrj3ijHExRi3PptPdQTXf3lVDY/O38Cp97zHr19aTmZ6B566fjw/P38Yb63YydP/3eJ1iGGlfoAuda+OfO3lmuGdlTspKq9W92ppU5OGdKdbpwRmLVA3axEvNNtvw8wuANYDDwB/A9aZ2XltHZj4HOtUTwCjM1PZuOeABupq4LWl23nm061847TjOP347lyek8lbPzidv149lljzJWhn/vl9Zn62NeISNeccP569hP2VNfzlqrFhMaVToHqnduB3l41i8dZi/vruOq/DabHSymr+/t46Trn7Xe6au5Lje3Zixk0TmPXNkznt+O58fWIWpw7uxp2vrGCtBl45ZPk2Xw+XURqgK2pE+zXD7LwCeqUkcepgzX0sbSc+NoYp2Rm8s2on+zSIpUjIBXJjw73AZOfcJOfc6cBk4L62DSuKnXceo37844BX75faD8PYWHRsLcgAy9XNGoBtxRX85PmljMlM5YdnH39oeWyMcdGYPrz2vVN56JocOiXF8aPnljDpj/P41yebqayu9TDqwD393y28s2oXPzl3KEN6Rd7IqheO7sNl4zL427trydu8z+twAlJcXsWf31rDKXe/yx/fWE123zSe+9bJPH3DBCYM7HpovfpW8k6JcXz3mUUR855qaxqgKypF7TXDrtJK5q3ZzWXjMoiNgttxJLxNzc2kutbxYn6h16GItDuBJMi7nHMNm3Q2ALvaKJ7oV1FB7MGDAa+eGJdIZkomG4pb3oJcf9G5tLC4xdtGm9o6x63P5lNTW8dfrhpLfCODXsTEGOeO7MXL35nIP6efQI+URH45Zxmn//E9/vfDjVRUhW9Ss25XGf/z6gpOHdyN6ScP8DqcY/abi0fQJ60Dtz6bT2kY93zYU3aQu19bxSl3v8sD76zlpOO68sp3J/LP68aT0z+90W16dE7iT9PGsGpHKXe/tirEEYenJQUaoCsKRe01wwsLC6nV3McSIkN7pTAqI1XdrEU8EEiCvNzM5prZdDP7GvAy8JmZXWZml7VxfIKvm/WxtCDXD9S1pED3If9j3jo+3biP314ykgHNDAJlZkwe2oPnv3Uy/3fDiWR168idr6xg4h/e5R/z1lN2sCZEUQemqqaOW59dRIf4WP40bUxEDzTWOSme+6/MprCogjteWuF1OF+wc38lv33Z9154eP56vjSsJ2/cehoPX5sb0EBTk4f24LpTBvDER5t4d9XOEEQc3jRAV1SKymsG5xyz8wrI6Z/OwO6dvA5H2olpuZms2L6f5dt0HScSSoEkyEnATuB0YBKwG+gCXARc2GaRySFZ6VnHdA8yQHa/NPK3Fgc3oAiTt7mI+95ey8VjfF14A2VmnDyoGzNuOolZ3zyJERmp/OH1VUz8w7t8vH5vG0bcMve/vYZlhfv5/WWj6ZmS5HU4rZY7oAvfnjyI5xYW8OqS7V6Hc8j//XcLp97zHk9+vIkLR/fhbf+96y3tzv7jc4f6BlGbtYRd+yvbKNrwV3Sgig17DpDdN83rUCS4ovKaYXFBCWt3lan1WELq4jF9SIiN0ZzIIiF2tGmeAHDOXReKQKRpA9MGUlhaSGVNJUlxLUuAxvVLZ+7SHezcXxkVyVNL7a+s5nszFtE7NYn/uXTkMU+xc8KALjx1/XgWby3mtlmLue6JT3nsqyd4Po3Sfzfs5R/vr+fK3L6cO7KXp7EE0y1nDGb+2j387IWljOufRu/UDp7G88R/NnLHyys4/fju/M+UkfTtknzM+0qKj+WvV4/lor99yA9nLebJ68ZHdKv/sVq0tQigyS7pEpmi9Zphdt5WkuJjuHB0b69DkXYkLTmBs4b35MX8bfz0vGFehyPSbgQyinWWmf3ZzJ43s5fqf0IRXFS68EL2nnRSizbJSvdN9bS5eHOLD1d/8blwc1GLt40Gv5qzjO0llfzlqrGkJLV+zsoxfdOYcdMEBnTtyPVPfsa81d7dWldSUc0PZi6mX5dkfnXRcM/iaAvxsTHcf2U2VTV13DZrMXV1zrNYHvtgA3e8vIJzRvTk0a/mtio5rje4Z2d+eeFwPli7h8f/0/LbJ6JB3uYiYmOMMZlpXociQRSN1wyV1bW8lL+N80b2pnMQziMiLTE1N5N9B6p4d1VU3MovEhEC6WI9B9gE/BXf6JT1P61mZuea2WozW2dmP2nkeTOzB/zPLzGzcYFuG7Zuu42tV17Zok1aM9XTiD6pJMTFkNcOE+TnFxYwJ38b3ztjcFBbqbp2SuSZGycwuEcnbnoqj7dXeHMv6a9fXMaO/ZXcf2U2HROb7QwScbK6deRXFw3nP+v2epZEPjhvHf/z6kouGNWbv315HAlxgfzLDMyXx/fj7OE9+cPrq9rlfOV5m4sY0SeFDgmRMx2ZBGQObXTN4JU3V+xkf2WNuleLJ04d1I0enROZnbfV61BE2o1ArvYqnXMPOOfec869X//T2gObWSzwd+A8YDhwtZkd2Qx2HjDY/3MT8I8WbBs1stJ8Lcgbi1ueJCTExTA6I5W8Le0rQd689wC/nLOM8f77WYMtvWMC/3fDBIb17sy3ns7j9WU7gn6Mo3kxv5A5+du45UuDGdsveruoXnVCX84a3pN7Xl/Nim2hna7sgXfWcs/rq7kkuw9/uSq70ZHPW8PM+MPlo+naMZFbnllEeVV4Df7Wlmpq61i8tYRxUfzebcdafM1QctCF9Ze4j3+4kZSkOBKD+AWZSKDiYmO4bFwm767axew1VWH7WcnbXMTf31sXtvGBYgyGcI8P/D3UOnVt1X2H5tzRuy6a2ZfxJahvAofmJ3LOLWzVgc1OAu5wzp3jf/xT/35/32Cdh4F5zrln/I9X4xv0Y0Bz2zYmNzfXLViwoDVht96kSRQXF5OWnx/wJs45kn+XTK9OvQ61JrfE3p2nUrxvLFlD/k5MTPhOVXRlxyt59sCzrd6PczEUbryS6qp0Mo/7F/HxpUGIrnG1tQls33IZByt60TPzVTqlrG2T4zSsm+qqzhRsuJb4xH1kDHgWM++6H4dCbU0SWzd8lZjYSjKz/o+YmM8TyWC9ZxpyDop2n0zRngl0Sl1Ojz5vtmkdVxzoy7bNU+mctowefd4K2n7bom6C5WBFDwo2XkOPjFfpnLr6mPYxIWMCd51xV4u3mzdvHpMmTTqmYzbGzPKcc7lB22GEO5ZrhsTeg13G9PsZ0Scl7Low76+oZtm2/RiQGB/D0zdMOGqPpGC/v7wSLeWA6CjLS/mF3DIjH4AYI+w+K6WV1Szftp86F3h8RUVFpKeH7kvSY4kxEMEsR1vFGKjmyuJ1fIGoj7Hwn7dycMfaYx7gJZB+maOAa4EvAXX+Zc7/uDUygIb9RQqAEwNYJyPAbQEws5vwtT7Ts2dP5s2b16qgWyu7uJja2toWx/HY6Mcoqyo7pmOuSU/iub1xTORqMjtWHdM+QqFLbBeu6HhFq/czryCFDZUpXHrcXoamnReEyI7u4DDHzDXVFBZcyIkD9zGia0XQj1FfN3UOntnanTjiuW6QkZ40LejHCkfrB1Ywc013uhZdx1n9Pu+OHKz3TD3n4P2CFDbsSWFMtzLOG5CKWRvXcUeYV1XKx9tHcWbXPgztEpz3T7DrJpgWlHWkAJjadRypiWOOaR+dyjsd0//zsrIyz88DUS6ga4aG5+aEXoOoc7B9bykHO4TXgHX7KnxfjjmgqrqOZ97+jNLjEppcP1reX9FSDoiOssxf//m1Wzh+VvZVOOqHCgk0vtraWvbsC10r5LHEGIhglqOtYgxUc2XxOr5AHIqxlWEFkiBfCgx0zgU7s2os9CObaZpaJ5BtfQudewR4BHwtyJ5/i5iWRnFxcUi/zdxdepDn7nqbvgNP46bTjgvZcVtq3rx5TJ00tVX7+Gj9Hu5e8F+uzO3LH6ZeEKTImndDVQ3XP/EZr2w0TssZE/R71err5h/z1rNlwSr+OHU003L7BvUY4S75peU88dEmbrngDE47vjsQnPdMPeccd726ko93bOSaCf347cUjQza69PW1dUx96GPeKYzjtmmnkZHW+lG7g1k3wbb8mUX0Kt7H7Vd//ZhHlj9W0dCaFOYCumZoeG5O7D3YJcXH8PB1R2+d9ULe5iK+8tgnVNfUER8Xw9VnnqAW5AgTDWXpnFXEK5s+oaq6joQw/Kwc+TkJJL5Qvy7HEmMgglmOtooxUM2Vxev4AlEfI811kW5GIAnyYiANCPbweQVAwyv8TGBbgOskBLCt+HXvnEj/rslhfb9AMBQdqOIHzy4mq2tHfn1xaG9JT06I45/Tx3PjUwu4ffZiauvquPKEfkE9xrLCEv781mrOH9WrXQ4W85PzhvKfdXv44azFvHHraXTp2HQrTks55/jNyyt44qNNTD95AL++aHhIE7f42BgeuCqb8//yAd+fkc8zN00gNoqnfsrbXERO//SQJ8cSEi2+ZkhPtGa7Lnslp386T98wgU827GXCwK5hGaNEv/r34TNvf9bslzReiITPiWJsvXCPDz6P8cT7i1uVFwaSIPcEVpnZZ3x+P5Fzzl3SmgMDnwGDzSwLKASuAr58xDovAd8xsxn4ulCXOOe2m9nuALaVBnL6pTN/7R6cc1F5Ueqc48fPLWHvgYM89rVTSE4I/ajOHRJieexruXzjX3n8+LmlVNc6rpnQPyj7PljruGXGIrp0TOCuKaOi8jVsTlJ8LPdflc2Uv/+Hnzy3hIevzQnKfuvqHL98cRlP/3cLN56axc/OH+ZJ/fbv2pE7p4zkBzMX8/f31nHLGYNDHkMo7CippLC4gusnZnkdirSNFl8zpCZaWF5o1cvpnx7W8Un7kNM/ndLjEsL2vRgJnxPF2HrhHh/4Yqwt29uq0XMDySJ+3eBvAyYCV7fmoADOuRoz+w7wBhALPO6cW25m3/Q//xAwFzgfWAeUA9cdbdvWxhQSV1zBrjVrSAvxYcf1T+f5RYVs3VdBv66tn8c13Pzfp1t4c8VOfn7+MEZmpHoWR1J8LI98NYdvP72QX8xZRk1tHdNPaX0i8OzqKjbsruHfXz+R9CC2nEaaEX1Suf2cIfxu7ipmLthKz1bur67O8dPnl/Lsgq18a9Jx/OicIZ5++XDp2AzeX7Obv7yzllMGdSWnfxfPYmkrC/0j6of7CVaOWZtcM4iIiIRKs3MW+KdnKAEuAJ4AzgAeCsbBnXNznXPHO+eOc87d5V/2kD85xvl82//8KOfcgqNtGxFuvpltU6aE/LD106nkbdkX8mO3tbU7S7nzlRWcOrgbXw+DVqnEuFge/EoO54zoyR0vr+CxD1o+f3VD767aybtbarhhYhYTB3cLUpSR64aJAzn5uK785uUV7DhQ1/wGTaitc9w2ezHPLtjKLWcM9jw5Bt/UT3dOGUmftCS+NyOf/ZXVnsbTFvI2F5EYF8Pw3ilehyJtoC2vGUREREKhyRZkMzseX9flq4G9wLP4poWaHKLYolN5OTGVlSE/7JBenemYEMvCzcVcOjZ67l+trK7lu88somNCHPdeMSZkgyo1JyEuhr99eRy3zsjnf15dSVVtHTdPano+5to6x/aSCgqKKti6r9z3u8j3e1lhCZmdjNvPHRLCEoSvmBjj3ivGcM5987nrvxW8sv2/9O3Sgcz0ZDLTfb/7dulA906JTSa8NbV1/GDmYl5avI0fnnU83w2j7swpSfH85aqxTHvoY37+wjIeuCrb88Q9mBZuKWJ0ZioJmlM2quiaQUREosXRulivAj4ALnLOrQMws++HJKpodv75jC4uhnPPDelhY2OMsf3So26grn99vJlVO0p5fHouPToneR3OYeJjY/jLVdnExRr3vL6ayuo6Thvc7QtJ8NaicrYXV1JT9/mAe2bQKyWJzPQOnDeyN+M77iUxLtbD0oSX3qkdeOxrJ/CnFz+l9GANby7fyd4Dhw+amxgXc1jCnJmeTF9/Ev3I/A28unQ7Pz53KN+aFH4ju4/rl86tZwzm3rfWcPUJfTl5UHT0HKisrmVZYYnuP45OumYQEZGocLQE+XJ83wa/Z2avAzNo9axS4qVx/dP527trKTtYQ6fE0A9i1RaeW1jA2H5pfGloa+9GbRtxsTH8+Yps4mJieOCdtTzwztpDz3XvnEjf9A6M7ZvOxWM+bwHtm55M77SkwxLiSJ/DsS2Mz+rCzdlJTJp0CgDlVTUUNmh5b/glxOKCYorLD++u/IsLhnHDqQO9CD0gN542kEfmb+C5hYVRkyAvKyyhutaR00/3H0chXTOIiEhUaDJLcs69ALxgZh2BKcD3gZ5m9g/gBefcm6EJUYIlp386dQ4Wby3mlCi44F61Yz+rdpTy20tGeB3KUcXGGH+cOppzR/YiPtYOJcJJ8WoRDqbkhDgG9+zM4J6dG32+tLKagiJfN/aUpDhOHNg1xBG2TFJ8LOeN6sWrS7bzP1NG0iEh8t8v9T1YxmmArqijawYREYkWgQzSdcA597Rz7kJ88w3nAz9p68Ak+LL7pmFG1HSznrNoG3ExxgWjensdSrNiYoyzhvdk0pAeDOrRScmxBzonxTOsdwpnDe8Z9slxvSljMzhQVcvbK3d6HUpQ5G0uYkDXZLp1SvQ6FGkjumYQEZFI16JRUpxz+5xzDzvnvtRWAUnbSe0Qz/E9OkdFglxX53gxv5DTju9OV11sS5SakNWVXilJzFlU6HUoreacY+GWIrUetyO6ZhARkUikYURDbfp0doR4gK6GxvVPY+GWIuoaDAgVif67cR/bSyqZMjbD61BE2kxMjHFJdh/eX7ObfUcMQhZptuwrZ09Z1aEp50RERETCkRLkUPM6Qe6XTmllDet3l3kWQzC8mF9Ix4RYzhoWnoNziQTLlLEZ1NQ5Xl2yzetQWmXhFl/PlRy1IIuIiEgYU4Icanv2EF9S4tnh6y9OI7mbdWV1La8u3c45I3tFxcBFIkczrHcKQ3p25oUI72adt7mITolxHN/EIGoiIiIi4UAJcqhNncqIX//as8NndetIenJ8RCfI763aRWllDZeqe7W0E1PGZrBwSzGb9x7wOpRjlre5mLH90oiN0cw/IiIiEr6UILczZkZO/3TytkRugjwnv5DunRM5+bjIn6pKJBAXZ/cB4MX8yOxmXVpZzeod+3X/sYiIiIQ9Jcjt0Lj+6WzYfSAiB/0pLq/ivVW7uXhMH7VESbuRkdaBE7O6MGdRIc5F3gB7i7eWUOd0/7GIiIiEPyXI7VCOvxVnUQS2Is9duoOq2jp1r5Z259KxGWzYc4Clhd6NYXCs8jYXYQbZ/dK8DkVERETkqJQgt0OjM9OIi7GIvA95Tn4hx3XvyIg+KV6HIhJS543qTUJsTEQO1pW3pYjje3QmJSne61BEREREjkoJcqh961sUXnyxpyF0SIhleJ+UQ9OuRIqConI+3biPS8dmYKbu1dK+pHaI50tDe/Dy4m3U1NZ5HU7A6uoci7YUMU7dq0VERCQCKEEOtSuvZPeXvuR1FIzrl87irSVUR9CFdv0ARZdkq3u1tE9Txmawp6yK/6zf63UoAVu3u4zSyhrdfywiIiIRQQlyqG3dSuKuXV5HQU7/dCqqa1m1vdTrUALinGPOokJy+6fTt0uy1+GIeGLy0O6kJMUxJ4K6WdffyqEEWURERCKBEuRQu/Zahv3ud15HcehiNW/zPo8jCcyK7ftZu6uMKRqcS9qxxLhYLhjdmzeW76C8qsbrcAKSt7mILh0TGNBVX2yJiIhI+PMkQTazLmb2lpmt9f/+QtOCmfU1s/fMbKWZLTez7zV47g4zKzSzfP/P+aEtQeTrk9aB3qlJ5G0p9jqUgMxZVEh8rHHBqN5ehyLiqSnZGZRX1fLWip1ehxKQhZuLGNcvXeMGiIiISETwqgX5J8A7zrnBwDv+x0eqAX7onBsGTAC+bWbDGzx/n3Mu2/8zt+1Djj7j+qezMAJGsq6tc7y0eBunH9+D9I4JXocj4qkTBnShT2pSRIxmve9AFRv2HFD3ahEREYkYXiXIlwBP+v9+Ephy5ArOue3OuYX+v0uBlYD61wbRuH7pFBZXsKOk0utQjuqTDXvZuf+g5j4WAWJijEvGZvDB2j3sKTvodThHVf8F3DjNfywiIiIRwpxzoT+oWbFzLq3B4yLnXJNNDGY2AJgPjHTO7TezO4DpwH5gAb6W5kabQs3sJuAmgJ49e+bMmDEjSKU4Ntm33kptbS1L//pXT+MA2FBcy28/qeTb2Ymc0CvO63AoKyujU6dOX1j+2NKDLNhRwwNfSiYhtn1202yqbtq79lovBaV1/OI/FXxlWAJn9W98buFwqJvZa6p4bWM1D56ZTGKYfHaDXS+TJ0/Oc87lBm2H7UTDc3P37t1zZs6c6XFEwREOn7tgiJZygMoSrqKlLNFSDoiusrT23NxmCbKZvQ30auSpnwNPBpogm1kn4H3gLufc8/5lPYE9gAPuBHo7565vLqbc3Fy3YMGClhYluF5+maVLlzLqZz/zNg6gqqaOUXe8wTUT+vPLC4c3v0EbmzdvHpMmTTpsWWV1Lbn/8zbnjezFH6eN8SawMNBY3Uj7rpfz/vIBCXExvPjtUxp9Phzq5sqHP6ayupYXvzPR0zgaCna9mJkS5FYaMmSIW716tddhBEU4fO6CIVrKASpLuIqWskRLOSC6ytLac3ObNRs6585s6jkz22lmvZ1z282sN9DovEdmFg88Bzxdnxz7972zwTqPAq8EL/I2dtFF7O3c2esoAEiIi2FMZtqhaVjC0dsrd1J2sEbdq0WOcOnYPvxu7io27jlAVreOXofzBdW1dSwuKObq8f28DkVEREQkYF7dg/wS8DX/318DXjxyBfMNefq/wErn3J+PeK7hUMaXAsvaKM7gW72aDlu2eB3FIeP6p7N8WwmV1bVeh9KoOYu20TMlkRMHdvU6FJGwcvGYDMwI2zmRV27fT2V1nQboEhERkYjiVYJ8N3CWma0FzvI/xsz6mFn9iNSnANcCX2pkOqd7zGypmS0BJgPfD3H8x+4b32DIn//c/HohktM/nepax9LCEq9D+YJ9B6qYt3oXl2RnEBsTHvcvioSLXqlJnDSwK3PyC/FiLInm1PdMUYIsIiIikcSTkZmcc3uBMxpZvg043//3h0CjWZFz7to2DbAdqR9dNm9zEScM6OJtMEd4del2auocU7LVvVqkMVPGZvCj2UvI31rM2H7hlYjmbS6iT2oSvVM7eB2KiIiISMC8akGWMNG1UyIDuiaH5X3ILy4q5PienRjWOzzu2RYJN+eO7EVCXExYdrNeuLmIsWo9FhERkQijBFkY1z+dRVuKwqqb5pa95SzYXMSUsRn4bkcXkSOlJMVz1rCevLJkO9W1dV6Hc8j2kgq2lVSSE2at2iIiIiLNUYIs5PRPZ09ZFVv2lXsdyiEv5vtaxC5R92qRo5oyNoO9B6r4cO0er0M5ZOHmYkD3H4uIiEjkUYIcar/4BZuvDa9bqOsvYsOlm7Vzjjn5hYzP6kJGmu5fFDma04/vTlpyPC+EUTfrvM1FJMXHMLxPitehiIiIiLSIEuRQO/NMinJyvI7iMIN7dKZzYlzYJMjLCvezfvcBzX0sEoCEuBguGNWbN1fsoOxgjdfhAJC3pYjRmWnEx+oUIyIiIpFFVy+hlp9Pp3XrvI7iMLExRna/tLBJkF9YVEhCbAznj+zd/MoiwqVjM6isruPN5Tu8DoXK6lqWF5aoe7WIiIhEJCXIoXbrrQz629+8juILcvqns3pnKaWV1Z7GUVNbx8tLtjF5aHdSk+M9jUUkUuT0TyczvUNYdLNeUlBCTZ3TAF0iIiISkZQgCwDj+qXjHORvLfY0jo/W72V36UF1rxZpATNjSnYG/1m3h12llZ7GUt8TZax/jnURERGRSKIEWQDI7peG2eejz3plzqJCOifFMWlID0/jEIk0U8b2oc7By4u3exrHwi1FZHXrSNdOiZ7GISIiInIslCAL4JtPdUjPzuRt8e4+5IM1jjeW7+CCUb1Jio/1LA6RSDSoR2dGZqQcmiLNC845Fm4uYpy6V4uIiEiEUoIsh4zrn86izUXU1TlPjr9oVy0HqmqZou7VIsdkSnYGSwpKWLerzJPjb95bzt4DVRqgS0RERCKWEuRQ+93v2HDDDV5H0aicfumUHqxhrUcX1x9tr6FPahLjB3Tx5Pgike7iMX2IMTxrRa6//1gJsoiIiEQqJcihdvLJ7B850usoGlV/UevFdE97yw6ybE8tF2dnEBNjIT++SDTokZLEKYO6MSe/EOdC3xMkb0sRnRPjGNyjU8iPLSIiIhIMSpBD7aOPSFm2zOsoGtW/azJdOyZ4kiC/mL+NOodGrxZppSnZGWzdV8GaorqQH3vh5iLG9k/Xl1wiIiISsZQgh9rPfsbAxx7zOopGmRlj+6WzMMQDdVVW1/Lw/PUMSothSK/OIT22SLQ5d2QvunRM4KX1VSE97v7KalbvLGWcpncSERGRCKYEWQ6T0z+djXsOsO9A6C6un/p4Ezv3H2Tq8QkhO6ZItOqYGMfNk45j+d46Plq3J2THXby1GOd0/7GIiIhENiXIcpj6i9uFIepmvb+ymgfnree047sztIumdhIJhmsm9KdLkvGHN1aH7F7kvM1FmEF237SQHE9ERESkLShBlsOMzkwlLsZCNh/yY/M3UFxezY/OGRKS44m0B0nxsVwyKJ7FW4t5c8XOkBwzb3MRQ3p2pnNSfEiOJyIiItIWPEmQzayLmb1lZmv9vxvtk2dmm8xsqZnlm9mClm4vLZcUH8uIjNSQDNS1p+wgj324kQtG9WZkRmqbH0+kPZnYJ46B3TvypzdWU9vGc5vX1jnytxSre7WIiIhEPK9akH8CvOOcGwy843/clMnOuWznXO4xbh9e7r+fdd/5jtdRHFVOv3QWby2murZtR8H9+3vrOFhTxw/OPr5NjyPSHsXGGD88awhrd5UxZ1Hbzou8dlcppQdrlCCLiIhIxPMqQb4EeNL/95PAlBBv753sbMoGDfI6iqPK6Z/OwZo6Vmzb32bHKCgq5+lPtjB1XCbHddecqSJt4byRvRiZkcJ9b6+hqqbtvvCq73GiBFlEREQinYVqAJfDDmpW7JxLa/C4yDn3hSsrM9sIFAEOeNg590hLtvc/dxNwE0DPnj1zZsyYEcyitFh6Xh4VFRVUTpzoaRxHs6+yjh/Mq+DLQxM4e0Db3E/4v0sP8vG2Gv5wWge6dvB9T1NWVkanTkqWG6O6aZzqpWn1dbN0dw335h3kmmEJnNm/bT7Pjy45yJI9NTwwORmz8J4DOdjvmcmTJ+cd0cNJAtDw3Ny9e/ecmTNnehxRcETL/6RoKQeoLOEqWsoSLeWA6CpLa8/NbZYgm9nbQK9Gnvo58GSACXIf59w2M+sBvAV81zk3vyUJckO5ubluwYIFza3WtiZNori4mLT8fG/jaMbJv3+H7H5pPPiVnKDve92uMs6+732mn5zFry4afmj5vHnzmDRpUtCPFw1UN41TvTStvm6cc1z1yCes332A+T+aRHJCXFCP45xj0p/mcXzPzjz61fDPE4P9njEzJcitNGTIELd69WqvwwiKaPmfFC3lAJUlXEVLWaKlHBBdZWntubnNulg75850zo1s5OdFYKeZ9Qbw/97VxD62+X/vAl4AxvufCmh7OXZnDu/Ja8t2MG918Kv2z2+tpkN8LN+efFzQ9y0ihzMzfnTuEPaUHeSf/9kU9P0/9P4GNu8t56xhPYO+bxEREZFQ8+oe5JeAr/n//hrw4pErmFlHM+tc/zdwNrAs0O2ldX5y3lCG9krhlmcWsXHPgaDtd2lBCXOX7uDrpw6ka6fEoO1XRJqW078LZwztwcPvr6ekvDpo+523ehf3vLGKC0f3ZlpuZtD2KyIiIuIVrxLku4GzzGwtcJb/MWbWx8zm+tfpCXxoZouBT4FXnXOvH217CZ7khDgeuTaH2BjjpqcWUHawJij7veeNVaQnx3PjqVlB2Z+IBOa2c4ZQerCGh+avD8r+Nu45wC3PLGJorxTumTo67O89FhEREQmEJwmyc26vc+4M59xg/+99/uXbnHPn+//e4Jwb4/8Z4Zy7q7ntJbj6dknmb18ex4Y9B/jBs/nUtXIu1Y/X7+WDtXu4edIgOie1zWBBItK4Yb1TuHhMH/75n43s2l/Zqn2VHazhpqcWEBtjPHJtTtDvaxYRERHxilctyO3Xww+z+gc/8DqKgJ0yqBs/O38Yb67Yyd/eW3fM+3HOcc8bq+iVksS1J/UPYoQiEqgfnHU8NbWOv7577J/lujrHD57NZ8OeA/zty+Po2yU5iBGKiIiIeEsJcqgNGUJFv35eR9Ei158ygMvGZvDnt9bw1oqdx7SPt1fuYtGWYm45YzBJ8bFBjlBEAtG/a0euPKEvz3y6hS17y49pH397bx1vrtjJz84fximDugU5QhERERFvKUEOtZdfputHH3kdRYuYGb+7bBSjMlL5/rP5rNtV1qLt6+ocf3pjNQO6JmsgHxGP3XLGYOJijfvfXtPibd9asZM/v7WGy8ZmcP0pA4IfnIiIiIjHlCCH2r330nfmTK+jaLGk+FgevjaHxLgYbnpqAfsrAx8J96XF21i9s5QfnD2E+Fi95US81DMlia+dPIAX8gtZvaM04O3W7Srj+8/mMyojld9dNkqDcomIiEhUUrYiAeuT1oEHvzKOLfvKuXVGYIN2VdXU8ee31jC8dwoXjuodgihFpDnfOv04OiXE8ac3Vwe0/v7Kam56agGJcTE8fG2ObpMQERGRqKUEWVrkxIFd+fVFw3l31S7uC6CL5rMLtrJlXzm3nzOEmBi1OImEg7TkBG46bSBvrdjJwi1FR123rs5x64x8tuwr58GvjKNPWocQRSkiIiISekqQpcWumdCfK3P78td31/Ha0u1NrldRVctf31nLCQPSmTSkewgjFJHmXD8xi26dEvjj66txruneIPe9vYZ3V+3i1xcN58SBXUMYoYiIiEjoKUGWFjMzfjtlBGP7pfHDWYubvI/xiY82sav0ILefM1T3K4qEmY6JcXx78iA+3rCXD9ftaXSd15Zu56/vruPK3L5cM0HTs4mIiEj0U4Icav/6Fyt/9jOvo2i1xLhYHromh46Jcdz41AKKy6sOe76kopqH3l/PpCHdGZ/VxaMoReRovnxiPzLSOvDHN77Yirx6Ryk/nLWYsf3S+O2UEfqSS0RERNoFJcih1rcvB3v08DqKoOiZksRD1+SwvaSC7z6ziNoGg3Y9Mn89JRXV3Hb2EA8jFJGjSYyL5dYzB7OkoITXl+04tLy4vIobn1pAx8Q4Hromh8Q4DcolIiIi7YMS5FB79lm6v/uu11EETU7/dO68ZCQfrN3DPW+sAmB36UEe/3ATF47uzciMVI8jFJGjuWxcJoN6dOJPb66mts5RW+f47jOL2F5SwUPX5NAzJcnrEEVERERCJs7rANqdf/yDjOJi+O1vvY4kaK4a349l20p4+P0NjOiTysLNRVTV1vFDtR6LhL3YGOO2s4/nm/9eyPMLC1i3u4wP1u7h7stGkdM/3evwREREREJKCbIExa8uHMHqHaX8aPZiauscV+RmktWto9dhiUgAzhnRi9GZqfz2lRWUVtZwzYR+XDW+n9dhiYiIiISculhLUCTExfDgV3JI65CAmXHLGYO9DklEAmRm3H7OEEorazhhQDq/unCE1yGJiIiIeEItyBI03TsnMvtbJ7G79CC9Uzt4HY6ItMDEQd146vrxjOmbRkKcvjsVERGR9kkJsgRVZnoymenJXochIi1kZpx2fHevwxARERHxlBLkUJs9m+X/+Q+neB2HiIiIiIiIHMaTfnRm1sXM3jKztf7fXxgq1cyGmFl+g5/9Znar/7k7zKywwXPnh7wQx6pbN6pTNfWRiIiIiIhIuPHqRrOfAO845wYD7/gfH8Y5t9o5l+2cywZygHLghQar3Ff/vHNubiiCDoonnqDX6697HYWIiIiIiIgcwasE+RLgSf/fTwJTmln/DGC9c25zWwYVEkqQRUREREREwpJXCXJP59x2AP/vHs2sfxXwzBHLvmNmS8zs8ca6aIuIiIiIiIi0hDnn2mbHZm8DvRp56ufAk865tAbrFjnnGk1yzSwB2AaMcM7t9C/rCewBHHAn0Ns5d30T298E3ATQs2fPnBkzZhxzmYIh+9Zbqa2tZelf/+ppHOGorKyMTp06eR1GWFLdNE710jTVTeOCXS+TJ0/Oc87lBm2H7UTDc3P37t1zZs6c6XFEwREtn7toKQeoLOEqWsoSLeWA6CpLa8/NbZYgH/WgZquBSc657WbWG5jnnBvSxLqXAN92zp3dxPMDgFeccyObO25ubq5bsGBBKyIPgkmTKC4uJi0/39s4wtC8efOYNGmS12GEJdVN41QvTVPdNC7Y9WJmSpBbaciQIW716tVehxEU0fK5i5ZygMoSrqKlLNFSDoiusrT23OxVF+uXgK/5//4a8OJR1r2aI7pX+5PqepcCy4IanYiIiIiIiLQ7Xs2DfDcw08y+DmwBpgGYWR/gMefc+f7HycBZwDeO2P4eM8vG18V6UyPPh6+5c1kyfz6neR2HiIiIiIiIHMaTBNk5txffyNRHLt8GnN/gcTnQtZH1rm3TANtScjJ1SUleRyEiIiIiIiJH8KqLdfv14IP0mTPH6yhERERERETkCF51sW6/Zs6kR3Gx11GIiIiIiIjIEdSCLCIiIiIiIoISZBERERERERFACbKIiIiIiIgIoARZREREREREBABzznkdQ8iY2W5gs9dxAN2APV4HEYZUL01T3TRO9dI01U3jgl0v/Z1z3YO4v3bHzEqB1V7HESTR8rmLlnKAyhKuoqUs0VIOiK6yDHHOdT7WjdvVKNbhchFjZgucc7lexxFuVC9NU900TvXSNNVN41QvYWl1tLwm0fL+ipZygMoSrqKlLNFSDoi+srRme3WxFhEREREREUEJsoiIiIiIiAigBNkrj3gdQJhSvTRNddM41UvTVDeNU72En2h6TaKlLNFSDlBZwlW0lCVaygEqyyHtapAuERERERERkaaoBVlEREREREQEJcgiIiIiIiIigBLkkDKzc81stZmtM7OfeB2Pl8zscTPbZWbLGizrYmZvmdla/+90L2P0gpn1NbP3zGylmS03s+/5l7frujGzJDP71MwW++vlN/7l7bpeGjKzWDNbZGav+B+3+7oxs01mttTM8uunfFC9hI9IPice5X/1HWZW6H/P5ZvZ+V7HGoho+ayY2ZAGdZ9vZvvN7NZIeF1ael1kZj/1f3ZWm9k53kTduCbK8kczW2VmS8zsBTNL8y8fYGYVDV6bhzwLvBFNlKXJ91MEvi7PNijHJjPL9y8P29flWK6VW/q66B7kEDGzWGANcBZQAHwGXO2cW+FpYB4xs9OAMuAp59xI/7J7gH3Oubv9F0vpzrkfexlnqJlZb6C3c26hmXUG8oApwHTacd2YmQEdnXNlZhYPfAh8D7iMdlwvDZnZD4BcIMU5d6E+T76LfiDXObenwbJ2Xy/hINLPiUf5X30FUOac+5OX8bVUNH5W/O+xQuBE4DrC/HVpyXWRmQ0HngHGA32At4HjnXO1HoV/mCbKcjbwrnOuxsz+AOAvywDglfr1wk0TZbmDRt5Pkfi6HPH8vUCJc+634fy6tPRa+VheF7Ugh854YJ1zboNzrgqYAVzicUyecc7NB/YdsfgS4En/30/ie7O3K8657c65hf6/S4GVQAbtvG6cT5n/Ybz/x9HO66WemWUCFwCPNVisummc6iU8RPQ58Sj/q6NJpH9WzgDWO+c2ex1IIFp4XXQJMMM5d9A5txFYh+8zFRYaK4tz7k3nXI3/4SdAZsgDOwZNvC5NibjXpZ6/IeIKfIlkWDuGa+UWvy5KkEMnA9ja4HEB0Xcyba2ezrnt4HvzAz08jsdT/m/vxgL/RXVT34U4H9gFvOWcU7187n7gR0Bdg2WqG9+XKG+aWZ6Z3eRfpnoJD1FzTjzifzXAd/zdSB+3COiW7BeNn5WrOPxiPxJfl6Zeg0j//FwPvNbgcZb5bhF638xO9SqoFmrs/RTJr8upwE7n3NoGy8L+dQnwWrnFr4sS5NCxRpapf7s0ysw6Ac8Btzrn9nsdTzhwztU657Lxfes83szCrtuPF8zsQmCXcy7P61jC0CnOuXHAecC3/d3LJDxExTmxkf/V/wCOA7KB7cC93kXXIlH1WTGzBOBiYJZ/UaS+Lk2J2M+Pmf0cqAGe9i/aDvRzzo0FfgD8n5mleBVfgJp6P0Xs6wJczeFfKIX969KCa+UWvy5KkEOnAOjb4HEmsM2jWMLVTv99BfX3F+zyOB5P+O+xfQ542jn3vH+x6sbPOVcMzAPORfUCcApwsf8ewhnAl8zs36hucM5t8//eBbyAr0tVu6+XMBHx58TG/lc753b6v8yrAx4ljLpXHk0UflbOAxY653ZC5L4uNP0aROTnx8y+BlwIfMX5B0Hyd3vd6/87D1gPHO9dlM07yvspUl+XOHxjujxbvyzcX5cWXiu3+HVRghw6nwGDzSzL/83mVcBLHscUbl4Cvub/+2vAix7G4gn/PSD/C6x0zv25wVPtum7MrLt9PuJlB+BMYBXtvF4AnHM/dc5lOucG4Pu/8q5z7hraed2YWUf/4B2YWUfgbGAZ7bxewkhEnxOb+l9df3Hmdym+91xYi9LPymGtYZH4uvg19Rq8BFxlZolmlgUMBj71IL6Amdm5wI+Bi51z5Q2Wd/cPqIaZDcRXlg3eRBmYo7yfIu518TsTWOWcK6hfEM6vyzFcK7f8dXHO6SdEP8D5+EbtXA/83Ot4PK6LZ/B136jG983O14GuwDvAWv/vLl7H6UG9TMTX7WMJkO//Ob+91w0wGljkr5dlwK/8y9t1vTRST5PwjTrZ7usGGAgs9v8sr/+f297rJZx+IvmceJT/1f8ClvqXv4RvpFXP422mLFH1WQGSgb1AaoNlYf+6tPS6CPi5/7OzGjjP6/gDKMs6fPeB1n9eHvKve7n/fbcYWAhc5HX8AZSlyfdTpL0u/uVPAN88Yt2wfV2O8v83aJ8XTfMkIiIiIiIigrpYi4iIiIiIiABKkEVEREREREQAJcgiIiIiIiIigBJkEREREREREUAJsoiIiIiIiAigBFlERERE2jkzqzWz/AY/A7yOKRjMbLqZ7Tazx/yPJ5mZM7OvN1hnrH/Zbf7HT5jZ1CP2U3aUY3Tw11mVmXVrq7KIhIoSZJEIZGZdG5zEd5hZof/vMjN7sA2O94SZbTSzbx7j9u/5Y8sNdmwiIiJBUOGcy27ws6n+CfOJ5GvmZ51zNzR4vBS4ssHjq/DNd3tMnHMVzrlsYNux7kMknETyh12k3XLO7a0/iQMPAff5H3dyzt3cRoe93Tn30LFs6JybDCwIcjwiIiJtwswGmNlK/5fOC4G+Zna7mX1mZkvM7DcN1v25ma02s7fN7JkGLbHz6r8YNrNuZrbJ/3esmf2xwb6+4V8+yb/NbDNbZWZPm5n5nzvBzD4ys8Vm9qmZdTazD8wsu0Ec/zGz0QEUbwuQZGY9/fs/F3gtwHr5bYMv6AvN7J+BbCcSSZQgi0QR/8n1Ff/fd5jZk2b2ppltMrPLzOweM1tqZq+bWbx/vRwze9/M8szsDTPrHcBxnjCzB/wn6w31XbHMrLeZzfefOJeZ2altW2IREZGg6NAg8XvBv2wI8JRzbqz/78HAeCAbyDGz08wsB18L7FjgMuCEAI71daDEOXeCf/0bzSzL/9xY4FZgODAQOMXMEoBnge8558YAZwIVwGPAdAAzOx5IdM4tCbC8s4FpwMn4vgA4eMTzf2zY5bx+oXPuV/4v508H9gJ/C/B4IhEjzusARKRNHQdMxnei/Ri43Dn3I//J/wIzexX4K3CJc263mV0J3AVcH8C+ewMTgaHAS/hOtl8G3nDO3WVmsUBy0EskIiISfPXdhAFfCzKw2Tn3iX/R2f6fRf7HnfAlzJ2BF5xz5f7tXgrgWGcDoxvc55vq31cV8KlzrsC/r3xgAFACbHfOfQbgnNvvf34W8Eszux3fefuJFpR3Jr6keyjwDL5EuaHbnXOz6x80vAfZ3+r8NL7ea3ktOKZIRFCCLBLdXnPOVZvZUiAWeN2/fCm+k+4QYCTwlr8XVyywPcB9z3HO1QErzKynf9lnwOP+1uk5zrn8oJRCREQk9A40+NuA3zvnHm64gpndCrgmtq/h896aSUfs67vOuTeO2NckDm/JrcV3rW6NHcM5V25mbwGXAFcAAY/z4ZzbYWbVwFnA9/hignw0dwAFzjl1r5aopC7WItHtIIA/ka12ztWfYOv4/KS7vMGgJKOcc2e3ZN9+5j/OfOA0oBD4l5l9NRiFEBER8dgbwPVm1gnAzDLMrAcwH7jUfCM5dwYuarDNJiDH//fUI/b1rQa3Oh1vZh2PcuxVQB8zO8G/fmczq2/kegx4APjMObevhWX6FfBj51xtoBuY2YX4kupbWngskYihFmSR9m010N3MTnLOfew/WR/vnFt+LDszs/5AoXPuUf/JfhzwVBDjFRERCTnn3JtmNgz42N/jqgy4xjm30MyeBfKBzcAHDTb7EzDTzK4F3m2w/DF8vbgW+rsr7wamHOXYVf5boP5qZh3w3X98JlDmnMszs/1Ai1tznXMftXQb4IdAH+BTfz285Jz71THsRyRsKUEWacf8J92pwANmlorvf8L9wDElyMAk4HZ/t60yQC3IIiIS9pxznY54vAnfLUgNl/0F+Esj296Fb/wOzOyOBstXAQ1Hlf6Ff3kd8DP/T0Pz/D/123+nwd+fAROOPLaZ9cHXI/TNpsp2RKyHHaPB8oZxT2/k+U7+35MDOY5IJFOCLBLhjjipzcN/4mu43P+4UxPb5OPrFt2SY05vbN/OuSeBJ1uyLxEREWk5/21MdwE/8CfdjakAzjOzx46YCzmYcXTANxBoPL5buEQimn1+S6KISOPM7C/4Rt38y7HMhWxm7+EbUfsi59ziYMcnIiIiIhIMSpBFRERERERE0CjWIiIiIiIiIoASZBERERERERFACbKIiIiIiIgIoARZREREREREBID/B4xIoqZX4CK8AAAAAElFTkSuQmCC\n",
|
|
"text/plain": [
|
|
"<Figure size 1152x288 with 2 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"# Plot timedomain and frequency domain of the bandpassed signal\n",
|
|
"fig, (ax2, ax1) = plt.subplots(1,2, figsize=(16,4),sharey=True)\n",
|
|
"## Spectrum\n",
|
|
"plot_spectrum(ax1, fft, fft_freqs)\n",
|
|
"ax1.set_xlabel(\"Frequency [MHz]\")\n",
|
|
"ax1.grid(True)\n",
|
|
"ax1.margins(0.1, 0.1)\n",
|
|
"ax1.get_legend().remove()\n",
|
|
"ax1.set_xlim(0, 200)\n",
|
|
"\n",
|
|
"## Signal\n",
|
|
"ax2.set_title(\"Signal\")\n",
|
|
"ax2.margins(0.1, 0.1)\n",
|
|
"ax2.set_xlabel('Time [ns]')\n",
|
|
"ax2.set_ylabel('Amplitude')\n",
|
|
"ax2.grid(True)\n",
|
|
"ax2.plot(imp_time/ns, imp, label='Bandpassed')\n",
|
|
"ax2.plot(imp_time/ns, orig_imp, label='Original', zorder=-1, color='g')\n",
|
|
"ax2.axvline(imp_time[imp_offset[0]]/ns, color='r', ls='--')\n",
|
|
"ax2.legend()\n",
|
|
"\n",
|
|
"fig.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 5,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"### correlation\n",
|
|
"def correlation_and_lag(sig1, sig2, mode=\"same\", 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=False, fix_positive=False, **corr_kwargs):\n",
|
|
" corr, lags = correlation_and_lag(sig1, sig2, **corr_kwargs)\n",
|
|
" \n",
|
|
" lag = lags[np.argmax(corr)]\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",
|
|
" print(\"negative lag\")\n",
|
|
" lag += len(sig2)\n",
|
|
" \n",
|
|
" return lag, (corr, lags)\n",
|
|
"\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",
|
|
" return 2*np.pi* f_beacon * calc_lag /samplerate \n",
|
|
" \n",
|
|
" "
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"### Correlate a sinusoid with the bandpassed delta peak"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 6,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Known offset [ns]: 6.0\n",
|
|
"Known phase [ns]: 1.884955592153876\n",
|
|
"Expected lag: 3.0\n",
|
|
"Lag time multiplier [ns]: 2.0\n",
|
|
"Calc Lag: 1\n",
|
|
"Calc Lag time [ns]: 2.0\n",
|
|
"Phase lag: 0.6283185307179586\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"###########################################\n",
|
|
"# Correlate a sinusoid\n",
|
|
"\n",
|
|
"## define sinusoid properties\n",
|
|
"\n",
|
|
"# in-band or out-of-band\n",
|
|
"if True:\n",
|
|
" f_sine = 50 # MHz\n",
|
|
" corr_sig_timelength = 1/f_sine # us\n",
|
|
"else:\n",
|
|
" f_sine = 20 # MHz\n",
|
|
" corr_sig_timelength = 1/f_sine # us\n",
|
|
" \n",
|
|
"## define general correlation signal properties\n",
|
|
"corr_sig_samplerate = samplerate # MHz\n",
|
|
"corr_sig_samplelength = corr_sig_timelength * corr_sig_samplerate\n",
|
|
"corr_time = np.arange(0, corr_sig_timelength, 1/corr_sig_samplerate)\n",
|
|
"\n",
|
|
"corr_sig = sin_delay(f_sine, corr_time, t_delay=0, phase=-1/2*np.pi)\n",
|
|
"\n",
|
|
"## define the signals\n",
|
|
"if not True:\n",
|
|
" sig1 = imp\n",
|
|
" sig1_time = imp_time\n",
|
|
" sig1_offset = imp_time[imp_offset[0]]\n",
|
|
"else:\n",
|
|
" sig1_offset = 6*ns # us\n",
|
|
" sig1 = sin_delay(f_sine, corr_time, t_delay=sig1_offset)\n",
|
|
" sig1_time = corr_time\n",
|
|
"\n",
|
|
"\n",
|
|
"## correlate it and find the best correlation\n",
|
|
"calc_lag, (corr, lags) = find_best_lag(sig1, corr_sig, mode='full', fix_one_short=False)\n",
|
|
"\n",
|
|
"calc_lag += 1\n",
|
|
"\n",
|
|
"#####\n",
|
|
"lag_times = lags / corr_sig_samplerate # us\n",
|
|
"calc_lag_time = (calc_lag) / corr_sig_samplerate # us\n",
|
|
"\n",
|
|
"\n",
|
|
"\n",
|
|
"#(calc_lag_time - 1/corr_sig_samplerate)\n",
|
|
"#((calc_lag_time*corr_sig_samplerate - 1) - (calc_lag_time - 1)) /samplerate\n",
|
|
"\n",
|
|
"\n",
|
|
"\n",
|
|
"\n",
|
|
"#calc_lag_time = ((len(corr)-1)/2 - corr.argmax())/corr_sig_samplerate\n",
|
|
"\n",
|
|
"calc_phase_lag = 2*np.pi*f_sine*calc_lag_time\n",
|
|
"\n",
|
|
"print(\"Known offset [ns]:\", sig1_offset/ns)\n",
|
|
"print(\"Known phase [ns]:\", 2*np.pi*f_sine*sig1_offset)\n",
|
|
"print(\"Expected lag:\", sig1_offset*corr_sig_samplerate)\n",
|
|
"print(\"Lag time multiplier [ns]:\", 1/corr_sig_samplerate/ns)\n",
|
|
"print(\"Calc Lag:\", calc_lag)\n",
|
|
"print(\"Calc Lag time [ns]:\", calc_lag_time/ns)\n",
|
|
"print(\"Phase lag:\", calc_phase_lag)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 7,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"2.0\n"
|
|
]
|
|
},
|
|
{
|
|
"name": "stderr",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"<ipython-input-7-42b7d03bc798>:26: UserWarning: Matplotlib is currently using module://ipykernel.pylab.backend_inline, which is a non-GUI backend, so cannot show the figure.\n",
|
|
" fig.show();\n"
|
|
]
|
|
},
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABhxElEQVR4nO3ddXgUVxfA4d/dKFHiBHcJEIIUdyhSpLiUCuWjQCkUaItXoECLFS9arEDxQtHiFGhxCO7uMeKevd8fG1LSCAkkmWxy3+eZJ2RmduZkSebszNw5R0gpURRFUXIfndYBKIqiKNpQCUBRFCWXUglAURQll1IJQFEUJZdSCUBRFCWXMtU6gPRwdnaWRYsW1WTf164ZvpYpo8nu/5VtAlEUxVicPn3aT0rp8t/5RpUAihYtyqlTpzTZd8OGhq8HD2qy+39lm0AURTEWQoh7yc1Xl4AURVFyKaM6A1CAr7/WOgJFUXIIlQCMTdOmWkegKEoOoS4BGRtvb8OkKIryhtQZgLEZPNjwVd0EVhTlDWl2BiCEsBRCnBBCnBNCXBJCjNUqFkVRlNxIyzOAKKCxlDJUCGEGHBFC7JRSHtMwJkVRlFxDswQgDXWoQ+O/NYufVG3qVMTE6QkKjcbURGAvJUIIrUNSchC9Xs/JrX8RERJG3S7NMTU30zokJZNpeg9ACGECnAZKAj9LKY8ns04foA9A4cKFszbAbEKvl+y4+ISfdl/nR58QAKYtPMbwFmWpWsRB4+iUnMB7zzGOrF5OVNgDw/e71lO1VTfqdGmOTqfGiuRUIjs0hBFC5AU2AQOllBdTWq9atWpSqyeBM0RkMPjfTPPqEjh7P5Bfj93llk8YRRyt6OcQyuM4N0b5OuMXGsXbHm4MbV6G0m62mRe3kmNdP3GR/UuWEPb8OsLEBo9672KV154zO9YSF+2PuVUBanf5kKot62gdqvIGhBCnpZTVkszPDgkAQAjxHRAmpZya0jpGnQD8b8Gy1hDy+M235VKO8Pc2s9Q7lPkHbxEaHUuHygUZ3LQUhRyt3nz7So734PJtdi1YQtBTbxAWFK/akub93sPK1vD7Exsdw/7lm7l0YBP6uGCs8pak4Yc9KVfHS9O4ldeT7RKAEMIFiJFSBgoh8gC7gUlSym0pvcZoE0DAHVjWCmIi4J0pYG6T4qqPgyL4/cxDztwPxM7SlNae+WlQxgWzF6fh547C1bngXAo+2spzbJn31y2W/XMXJPSoWZjPGpXE2cYia342xaj43n/KzrlL8L1zDNBRwKMhLT79iLyuyV9KjAyLYPeCNdw8uQOpj8DOzZNmn/SiSMWSWRu48kayYwLwBJYDJhiGo66TUn6f2muMMgE8v2c4+EeHwkdbIV/FZFd7FBjBzL3X2XD6IVbmpvSpX5xedYthY/Gf2zQNG4JDMFR9BE6l4KMtYOXIk6AIZu69wbpTD8hjZkLvesXpXa8YtpbqRp4Cwf5B7Px5OQ8v7Qf0uBStSYv+vXAtki/Nr/9z7q88uLQfZCzORWrQol8v3Irnz9zAlQyR7RLA6zC6BBB4H5a2gqhgw4HavVKSVQLCopl74Ca/HrsHEj6oVYTPGpXE0do8+W2+qAb6yzewuju4lIEP/wArRwBu+oQybc81dlx4iqO1OZ81KkmPGoWxNDPJpB9Syc4iQ8PZtWA1t07tROojsXerRLO+/6Nw+eKvtT2/hz7snLsEn1v/ADryl61Py/49yevmmLGBKxlKJYCsFvjA8Mk/MtBwgM5fOdHisKhYFh+5w8JDtwmPjqVjlYIMfrs0BfLmSX27L5eDvrEX1nQHVw/4cDPk+fc0/tyDQKbsusaRm34UyJuHwU1L0aFKQUx0auhobhAbHcO+pb9z+a/N6ONCsHIoRaOPPqZsLc8M2f7Dq3fZvWAJzx+fAWFBsSrNaNHvfazsrDNk+0rGUgkgKwU9NBz8w5/Dh5ugQNWERVGxcaw+fp85B27iFxpN8/JufNWsDKXSOornv/0Aru+CNT0Ml5Y+2AR58iZa/cgNPybvusr5h0GUcrXhq+ZlaObhpp4hyKH0ej1H1vzJmZ3xo3isC1Kny4dUaVE7U/Z38+Rl9i1dQqj/VYSJNWXrtKVp786YW6RwBqtoQiWArBL8GJa+A+H+hgNywX/f8/1Xn/HtH5d4+DyCmsUdGd6iLJULp3Mcf3INYa7thLUfGC4xfbAJLO0SvURKyZ8XnzJl9zVu+4bhVSgvkzp6UiafGjqak1w7doHd82cTHfEYUwtnqrbuTu1Ob2fJOP7z+05w6LdlRIXex8TMgdpd/0f1Ng0zfb9K2qgEkBWCnxg++Yf6GA7Ehd5KWHTtaQht5xyhqJM1o1uVo14p59f7FP6iEqiXV+L5V7fDug8hfxV4f2OSJAAQG6dn45mHTNl1DUszE3YMqoedukmcIwQ+C2DpkM+QSDybdKbhR+9iapq1z3nq9XqO/3GAE7+vIDY6kI6jplC0UqksjUFJXkoJQD3il1FCnsLyNhD6zHAAfungHxkTx+erz2JracrK3jWoX9rl9S/BeHklPfgDlG0FnZfB4zOwqhNEhSRZxdRER9e3CrPgg6o8Dozg280pPnOnGBG9Xs/6cZPRx4XRsv8Imv6vY5Yf/AF0Oh212jeh29gfEDpTtkyfTHRkVJbHoaSdSgAZIdTHcPAPfgw9NkDhGokWT9x5lWvPQpjauRIutm84Pn/vXsOUnHJtoNMSeHgKVnWGqNBkV6taxJFBTUqz2fsxm84+fLN4FM3t+WUDwb4XKV6tLeXqemkdDm7F81O9/SfERDxh448/ax2OkgqVAN5UqK/h4B/0EHqshyK1Ei3ef/UZy/65S686xWhYxvXN9zd+vGFKice70PEXeHACfusC0WHJrvZZoxK8VdSBbzZf4r5/+JvHpWji7rkbXNz3G3nsivPukI+1DidB3S7NcS5ai8dX93Nq2yGtw1FSoBLAmwjzMxz8n9+D99ZB0cT1UnxCIhm6/jzl3O0Y3rJM1sVVoQN0WAj3j8JvXSE66QHe1ETH9K5eCAGfrzlLTJw+6+JTMkRUeBRbpk1C6MzoOHIEOtPs9axH56+HYGruzKHf5hLwxE/rcJRkqATwusL8YXlbeH4H3lsLxeolWqzXS75cd46w6FhmdfPCIqv/OCt2gvYL4d7fsLqboQzFfxR0sOLHDhXxfhDIrH03sjY+5Y1tnDiHmMinVG/fO1s+kWtla0XLgV8h4yJYP24Ser36kJHdqATwOsID4Nd3IeAWdF8DxRskWWXJ33c4fMOPb1p7pH2Mf0bz7Azt5sOdQ4anhpNJAq0989O5akHmHLjJsdv+GgSpvI6TWw/x5NoBXIrWpm6X5lqHk6LS1StQulYHQv2v8Oe81VqHo/yHSgBp1LBh/BD8Fwd/v+vQ7Tco0SjJuhcfBTHpz6s083Djveoa9zCo1BXazYXbBw0PjMVEJlllTNvyFHWyZshab4LCY7I+RiVdAh77cnj1z5haONPp68Fah/NKrT7/ACuH0lw5tJ6bp69oHY7yEpUA0sHG9DmsaA++Vw0H/5JNkqwTHh3LoDVncbQ2Z1JHz4x/4nbBAsOUHl7vQdvZcGsfrH0fYhMPzbO2MGVmNy98Q6IYuek8xvRsSG6j1+tZN24iMi6SlgOGJpRvzs50Oh2dvx6B0JmzfeYUIsOSnokq2lAJII1sTAOZ6tkefC5D11VQqmmy643bdoXbfmFM7+KFQ0oF3d5EmTKGKb2qfABtZsHNPYanhv+TBDwL5uWr5mXYceEp6049yKBglYy2c+5vhAVco3TtDpSuXl7rcNLMuaArdbp9SmyUDxsmzNI6HCWeSgBpISVfl+tNCZuL0GUFlG6W7Gp/XnzC6hP36degBLVLOmdOLFu3GqbXUfUjaD0DbuyCP0ckWdynXnFql3BizJbL3PJN/hkCRTs3T13h6uH1WDmUptXAD7QOJ91qvNuIfCXr8+zWYY7+nsKzLEqWUgkgLS7/QU2nPSy4/T2UaZHsKk+CIhi+8QKeBe0Z0rR05sXy00+G6XVV+xhqfganlhoeGHuJTieY1sULCzMdn68+S1Rs3BsGq2SUyNBwts+ajNBZ0PnrEUbbp7fjqIGYWrpxdP1CfO8/1TqcXM84f4uyUlQI/DmSGyEV2fSoT7KrxOklQ9Z6ExOnZ2a3ypibZvO3tdFIsM0H24ZAXGyiRfnsLZnU0ZNLj4P5afd1jQJU/mv9hFnERvlSp9unOBfMgAcKNWJpnYc2g4Yh9dFsmDARvfqQoalsfqTKBg5OhJAnTL8xnTiZfH2V+X/d4tjtAMa2LU8xZyOoh25hCy1+hKfn4eQvSRY3L5+PHjUKs/DQbQ7f8NUgQOVlRzfuxef2EfKVrE+Nd5OOOjM2xauUwaNBF8IDb7Jt9q9ah5OrqQSQmqcX4dg8qPoRl4PfSnaVs/efM23PdVp7utOpasEsDvANeLSDEk1g/3hDIbv/+LqVByVdbfhi3Tn8Q1VBL6343HvK0Q0LMLN0o9PoQVqHk2Ga9+uGjVM5bhzbxNWj57UOJ9dSCSAlej1s/8LQYKXJd8muEhoVy6A13uSzs2RC+4rG1WRFCEOD+rho2DUqyeI85ibM6laZoPAYhm1QQ0O1oI+NY+OEH5H6GFoPGoaF1RsWEsxGdDodnb8ZjjCx4s+5PxEenHzNKiVzqQSQEu9V8OA4vD0uod/uf337x0UePg9nRjcv7PNkUV39FSsMU0ZwKgH1voCLG+HWgSSLPfLbMaJlWfZd9WHFsXsZs08lzbbOXE540C08GnSleJUsrCWVRRzdnan//mfERfuzfvx0rcPJlVQCSE54AOz5FgrXgkrdk13lD+9H/H7mEQMbl+KtolnYELtQIcOUUeoMBsfisOOrJM8GAHxcpygNy7gwfvsVrj1N2mNAyRxX/znPzRObsXH2oHm/rlqHk2mqvVOPAuWa4HfvGEfW/Kl1OLmOSgDJ2fsdRAZBq2mQzHC7BwHhfL3pIlWLODCwccmsjW3tWsOUUcws4Z2p4H8T/k76gI4QgimdKmFnacrnq88SGaNGbWS28OAw/pz3E8LEiq7fDDfaIZ9p1WFEf8zz5OfEH7/w9PYjrcPJVXL2b9bruH8czvwKtfqDm0eSxbFxegav9QZgRlcvTE2y+C2cN88wZaSSTQw3hQ9PhYA7SRa72FowtXMlrj0LYeLOqxm7byWJdeOmERftT4P3B5A3n5PW4WQ6c0sL2n45HCnj2PjjRGJjY1/9IiVDqATwsrhYw41fuwLQIOmTsgCz99/k9L3nTOhQkUKO2b8OS5q1+BF0prBzGCRzw7dhGVd61SnGsn/usv/qMw0CzB0Ord6B//3jFCjXlKrv1NU6nCxTpGIJPJv2IDL4DlumLdE6nFxDJYCXnVgIzy5Ci4lgYZNkcUhkDLP336BjlYK0rZT96q+/Ebv80GgU3NgNV5IvNTG8ZRnKudvx1frz+AQnrSqqvJknNx9y6o8lmFvlp8PI/lqHk+Wa9OqAvZsnd05v49KhM1qHkyuoBPBC8GM4MAFKvm3orfsff+yMwbLdYQo5WjH2XeMpwpUu1fuCWwVDnaBk+glbmJowq5sX4dGxfLn+HHq9GhqaUWJjY/l94kQkcbz71QjMLTKhkGA2p9Pp6PLtMHSmNuxZOJ3Q52rQQWZTCeCFXaNAH2sYG5/MeP6xWy/xLDiSmd0qY2OR/BPBRs/E1HDjO/gR/DUx2VVKudnyTWsPDt/wU0NDM9D2mb8SGXKXSm+/T+HyxbUORzN2znlp/PEg4mIC2fDDDK3DyfFUAgC4uQ8ubYJ6X4JjsSSLT98L4Pczj+jboDhehfJmfXwv27DBMGWWwjWgyodwdC48u5zsKu9VL0zdks78tPuaeko4A/jce8rNk1uxdixD44/bax2O5io1rYF76Qb43z/Olb+9tQ4nR1MJICbSMAbesQTUSfqovV4vGbPlMvnsLPmsURYP+UyOs7NhykxNx4KlveGGeDJ9XIUQfNfGg7DoOKbtUQXj3tT2mfNBxtGyf/8cP+QzrVoP/gShs2Lf4gWql3AmUr9tf8+AgNvQ6icwTfqo/YbTD7nwKIgRLctiZZ4NLv0sW2aYMpOVI7w9Fu4fhXPJ93Et5WbLBzWLsPrEfS4/Ds7ceHKwS4fOEPDoFPnLNqRIxRJah5Nt2DnZU75RJ6LCHnBwxWv2v1BeKXcnAP9bcHgaVOiYbG/f4MgYJu+6StUiDrzrlU1G/WRFAgDweh8K1YA93xiejE7GkKalsc9jxpitl1StoNegj41j/9IFCJ01bQb11jqcbOftXh0xy+POuV1rVK2gTJJ7E4CUsGMomJhDswnJrjJ73w38w6IZ06a8cRV6ywg6neGGcEQg7Bub7Cr2VmZ82awMJ+4EsP3Ck6yNLwfY/+tmosMfUbFJZ2wc7bQOJ9vRmZrQ8MM+6ONC2Dpjsdbh5Ei5NwFc/sPQJL3x12DnnmTxLd9Qlv59ly5VC1GxoL0GAWYD+SpAzU/h9HJ4cDLZVbpXL0w5dzt+2H6FiGhVJiKtQp+HcH7vOsytCtCkp7rxmxLPxm+R170yDy/t48GVu1qHk+PkzgQQ3+WLfBXhreRPvcdvu0weMxO+ap7zqjCmS8MRYOsO25N2DwMw0QnGtPHgcVAkCw7d0iBA47Rt5mJkXBiNevZBZ2qidTjZWuvPPwWhY+fPGVwCRcmlCSC+yxetZxjGvv/Hgas+HLjmy+dNSuFim3NqsL+WhO5hF5LtHgZQo7gTrTzdmf/XLR4FRmRxgMbn/qXbPLpyAIcCVanQoKrW4WR7bsXzU7RyC0J8L3F6xxGtw8lRcl8CeKnLFwWrJVkcHatn3LbLFHex5qPaRbM+vlfZscMwZSWPd//tHhac/LX+kS3LIiX8sONK1sZmhHb+PA+Ejtaf99M6FKPRasCH6EzzcmTNYmKjY7QOJ8fIXQkgDV2+lv1zh9t+YXzT2iN7Nne3sjJMWenl7mG7Rye7SkEHK/o1KMH28084dts/a+MzIie2HiTU/wrFqrTEtWjSe09K8iyt8/BWm/eJjfLlzwVrtA4nx8iGR7hM5L0y1S5fPiGRzNp3k8ZlXWlUxlWDANNg7lzDlNUSdQ/bn+wq/RqUIL+9JWO3XiZO1QlKIjoqmqPrlmFi5kCrAR9pHY7Rqd2lGXnsi3P97z8IfKo+ZGQEzRKAEKKQEOKAEOKKEOKSECJzO16H+b+yy9fUXdeIio3j61blMjWUN7JunWHSwovuYduT7x6Wx9yEUa3KceVJMGtO3s/6+LK5XfNWExvtx1vvfpCj+vtmFZ1OR7O+nyJlFFtmLNA6nBxByzOAWOBLKWU5oCbwmRAiaQeWjLL3O4gMTrHL1/mHgaw//ZCP6xSjuEvSUtAK/3YPC7gFf89MdpVWFd2pXsyRqbuuERSurtW+EPDYl+vHtmCVtyS1OjbVOhyjVbJqOVxL1MH3zlGun7ikdThGT7MEIKV8IqU8E//vEOAKUCBTdnb/OJxdkWKXLyklY7ZcwsnaIsUWjw0bGqZcr2QTKN8eDk01lND4jxd1goIiYpi+V9UJemHrjAUgo2ne71NV7+cNtRnUB6GzZM/CeapO0BvKFr+JQoiiQGXgeDLL+gghTgkhTvn6+r7eDs6uSLXL12bvR5y5H8iwFmWwtTR7vX3kJs1/ABMz2JF897Dy+e3pVr0wK47d4/ozVdP96tHz+N07hlvJuhSvnMufK8kAed0cKVu3HZEhd1Uj+TekeQIQQtgAG4HBUsokVcWklAullNWklNVcXFxebydtZsHHO5Pt8hUWFcvEnVepVNCeTlUKvt72cxu7/NBoNNzck2L3sK+alcHa3ITvt17O1XWC9Ho9e3+Zj9Dloe2QvlqHk2M069sVU0s3Tm9fRWRouNbhGC1NE4AQwgzDwX+VlPL3TNuRTgcORZJdNPfgTZ4FR/Ftm/LodEZQ7+fgQcOktep9wK1iit3DHK3NGfJ2aY7c9GPP5dzbQ/jQb9uJCr2PR/0O2Dnn1TqcHMPU1JT67/0PfWwQ22Yv1zoco6XlKCABLAauSCmnaRHDff9wFh2+Q/vKBahaxEGLEIyXiSm0Tr172Ps1i1DK1Ybx268QGZP76gSFh4RzdsdqzCzz0fSTzlqHk+NUbl4bO9eK3PPexZObD7UOxyhpeQZQB/gAaCyE8I6f3snKAMZvv4ypTjCiZdms3O2bmTrVMGUHhaq/1D0s6YgMMxMd37bx4H5AOIuP3NEgQG1tm7kEfVww9d/vjalpNuglkQO9M/BTALbPUnWCXoeWo4COSCmFlNJTSukVP2VZjYMjN/zYffkZnzUqiZudZVbt9s1t22aYsosX3cO2Jd89rF4pF972cOPnAzd5GhSpQYDaeHT9Pg8u7MXerRJeb9fUOpwcq0DpwhSq2JSgZ+fw3nNM63CMjuY3gbUQE6dn7NZLFHa04n91k/YAVtLByhHe/h4eHINzvyW7ytetyhEbJ5n059UsDk47O2YbntZu9fmnGkeS87Ue1AudiR2HVv5CbGzSirVKynJlAlh17B43fEIZ3aoclmaqFO8b8+ph6B62O/nuYUWcrOldrxibzj7i9L3nGgSYtc78+Q/BPhcpUqkZ7iXVyLLMZmVrReV3uhMT+ZS9v2zQOhyjkusSQEBYNNP2XKduSWeaebhpHU7O8KJ7WGQQ7B2T7CqfNSqJq60F32+9hD4H1wmKjY7h8OrF6Eztaf15T63DyTXqv9cKC5siXP5rI8F+gVqHYzRyXQL4afc1wqLj+K6Nh3G2ecyTxzBlNy+6h51JvnuYtYUpI1qW5dzDIDacybkjNnYvWkds5DOqtuqBpU0WV23NxXQ6HU1790XqI9g6faHW4RiNXJUALj8OZvWJ+3xQswil3Gy1Duf17NxpmLKjhiPANj9sS757WDuvAlQunJfJf14jJDLn1QkKfBbA1SObsLQtSt1uLbQOJ9cpW8sT5yI1eXrzCLfPXtM6HKOQaxKAlJIxWy9hn8eMIU1Lax1OzvSie9izC3ByUZLFOp1gTJvy+IVGMXv/TQ0CzFxbZyxE6qN4u4+q96OVNoP7gjBj13xVJygtcs1v6Y4LTzlxJ4Avm5XB3sqI6/2MG2eYsiuPd6FkU9g/IdnuYZUK5aVz1YIs/fsOt32TPkFsrG6euoLP7b9xKVqL0tXLax1OruWY34XSNdsSHniTY7/v0zqcbC9XJICI6Dh+2HGFcu52dK9eWOtw3sy+fYYpuxICWk42dA/bNSrZVYa2KIOFqQnjt+eM9pF6vZ7dC+chhIWq95MNNP+0O6bmzpzY/CtR4Un7Vij/yhUJYMEhQ7PyMW08MHnNej/ZpQSPUXAqAfW+hEu/J9s9zNXWkoGNS7L/qg+7Lz3VIMCM9c+63UQE3aZ0nXfJm89J63ByPXMLc2p16UlczHO2zlyidTjZWq5IAGXz2fG/usWoUVz9cWaZOoP+7R4Wk/QJ4I/rFKOcux0jfr+AT7BxPyGc190F+3xetOjbTetQlHjV2zTE3q0S97x3cunQGa3DybZyRQJoUSEf37TOvGZjSjJe7h72z6wki81Ndczq5kV4dCxfrj9n1M8GVGhQld4zx2NqbsT3lnKgLt8ORWdqw56F0wl9rvpSJCdXJIAcxcnJMBmDV3QPK+VmyzetPTh8w48lf+e+YnFK5rJzzkvjjwcRFxPIunHZpIBiNqMSgLHZuNEwGYvmP4KJOewYmmz3sPeqF6aZhxuT/rzKxUdBGgSo5GSVmtagcMXmPH90mgO//qF1ONmOSgBK5rJzh0aj4OZeuLIlyWIhBJM6euJobc7na84SHq2KeSkZq/2wvlhYF+LMjuU8unZP63CyFZUAjM3IkYbJmLzoHrZzBEQlvRbrYG3O9C5e3PELY9y2yxoEqORkpuZmtB9u6Ae+adJEYqNz3lPor0slAGNz9KhhMiYvuoeFPIaDyXcPq13SmX4NSrD6xAN2Xkj6AJmivIkCZYpQtVVPosIesGnyAq3DyTZUAlCyRqHqUOUjODYv2e5hAF+8XZpKBe0Z8fsFHgdGZHGASk7X8IO2OBSoyv0Luzi397jW4WQLKgEoWafpGMiTN8XuYWYmOmZ2q0xMnJ4ha72JM+KhoUr21OWbrzAxy8v+pTNV2WhUAlCyUhq6hxV1tub7dytw/E4A8/+6lcUBKjmdjYMtb/cZgj42lLXfT871BeNUAjA2BQsaJmNV6T0oVDPF7mEAHasUoE2l/Ezbc52z93N+BzEla5WvX4ViVVsT/Ow8+5b8rnU4mlIJwNisXGmYjJVOZ7ghnEr3MCEE49tVIJ+dJYPWeBtH7wB9nNYRKOnQ9oteWNoV4/zeVdy7kHvPNFUCULKeW/l/u4clUywOwD6PGTO7efHweTjf/ZH8TeNs4+ZeWNAAQoy/sF1uYWpqSseRIxDChC0/TSI6MndWDVUJwNgMHmyYjF3DkeBSDtb0gLtHkl2lWlFHBjYuxe9nH7H57KMsDjCNbh0w/AwCwxPPitHIV7wA1d/tTXTEY36fOFfrcDShEkAaNWxomDTn7W2YjJ2FDXy0BewLwaoucO+fZFcb2LgkVYs48PXmi9z3D8/iIF/h9kFY3Q2cSsKHWww3uRWjUrdbC5yL1OTRlX2c2nFY63Cy3CsTgBDCSgjxjRBiUfz3pYQQrTM/NCXHs3GFj7aCXX5Y1RnuH0uyiqmJjhldvRDAoLVniY3LJqM27hyC37oZSl5/+Ic6+Buxzl8PwcTciUMrfybgiZ/W4WSptJwBLAWigFrx3z8ExmdaREruYusGPbeBbT5Y2REenEiySiFHKyZ0qMjZ+4HM2ndDgyD/4+4R+K0rOBQxfPK3dtY6IuUNWNlZ06L/l8i4CNaPm5SrhoamJQGUkFJOBmIApJQRGK54KkrGsM1nOBOwcYUVHeDhqSSrtK2Un45VCjLnwE2O3/bXIMh49/4xXLKyLxQfs4t2sSgZpmwtT0rVbEeo/xV2zV+jdThZxjQN60QLIfIAEkAIUQLDGYGihdKltY4gc9jlh4+2wbJWsKI9fLgZClRNtMrYd8tz6l4AQ9Z6s3NQfeytsrgBy/1jhktVdvn/TVgKADExMTx8+JDISOPt7layeQ3yvVUaqY/lwrnzRtngx9LSkoIFC2JmlrbYhUymRnuiFYR4G/ga8AB2A3WAnlLKg28WavpVq1ZNnjqV9NNhVnhxA1j1Bc5kQQ9h6TsQGWi4tp6/cqLF3g8C6TTvH5qVd+Pn96ogRBadjD44aUhMtm7Qc7vhrEVJcOfOHWxtbXFycsq6/5NMEBMVg/+j+whhgkvhwuhMjGecjJQSf39/QkJCKFasWKJlQojTUspq/33NK386KeUeoAPQE1gNVNPi4K/kEvYFDfcELO3h13bw5FyixV6F8vJFs9LsuPCU9aceZk1MD0/Dyg6Gyz0fbVUH/2RERkYa/cEfwMzCDBsHF6Q+hudPn2kdTroIIXByckrXWViKCUAIUeXFBBQBngCPgcLx8xQt9OljmHKyvIUNl4MsbOHXd+HphUSL+9YvQa3iTny35RK3fUMzN5ZHpw2f/K2cDDHZ5c/c/RkxYz/4v2DjYIeZhQ0xkaGEPg/WOpx0Se//QWpnAD/FTz8Dx4GFwKL4fyft8q1kjevXDVNO51DEcCZgZg3L28LTiwmLTHSCaV0rYWGmY9Aab8KiMqmL2OOzhoN/nryGWOwLJLva87BoVh5TnaZyEgd3N4TOjNDnvkRHRmsaiz5OT+Azf/SZMAQ6xQQgpWwkpWwE3AOqSCmrSSmrApWBmxkeiaL8l0NR6LkVzPLAr23h2b/dwtzt8zCpoycXHgXRYMoBlv9zl+jYDPwDeXLOcAnK0j7+4J+0AF9YVCyz992g/uQDfPvHRW76JO12pmStCRMmUL58eTw9PfHy8uL48eP07t2by5fT12lOZ6LD3tUNpCTg0X0CHj8jLjbxBw0bG5uMDD0JqdcT5BuA7727RIYGEBaY8WcjaRkFVFZKmXAOLqW8KITwyvBIFCU5jsUN192XtYLlbQwHY9dyADQvn4/f+9dm0s6rfLflEr8cuc2Xb5ehbaX86HRvcDni6QXDpScLW8Nln7yFEy2OjtWz+sR9Zu+/gV9oNM083PiqeRlKutq+yU+qvKGjR4+ybds2zpw5g4WFBX5+fkRHR/PLL7+81vYsrfPgkL8Qwb5+REcE43s/FEsrO+xcHNGZmGRw9P+SUhIaEEx4UABSxiJ05tg6uWJtn/EJJy0J4IoQ4hdgJYahoO8DVzI8EkVJiVOJf4eILm9jGIXjUgaAKoUdWNOnJodu+DFp51UGr/Vm/l+3GNaiDI3KuKb/uvTTi4ZLTmbWhsTjUCRhkV4v+ePcI6btuc6DgAhqFHNk4YdlqVLYISN/WqM3duslLj/O2E+rHvnt+K5N+VTXefLkCc7OzlhYWADg7Gx4QK9hw4ZMnTqVatWqsXjxYiZNmkT+/PkpVaoUFhYWzJkzh549e2JnZ8epU6d4+vQpkydPplOnTsTExdDt4w/x9/cnKjKS4YMH0bJZcyxt82boz/dCWFAIoQH+SH0MQphi7eCKjYNdpt1fScsYp4+BS8AgYDBwOX6eogUvL8OU2ziXNByQhc6QBPz+fSJYCEGD0i5sG1iXWd0rExETR69lp+iy4Cgn7ybfcyBZzy4bLjWZ5TFcenI0DKWTUrL/6jPemXWYIWvPYWdpxvJe1VnTp6Y6+GcjzZo148GDB5QuXZr+/fvz119/JVr++PFjxo0bx7Fjx9izZw9Xr15NtPzJkyccOXKEbdu2MWKEoYm8paUlmzZtwtvbm8NHjjBu8lSk0BER7I+UkhD/IF41lD4tIkLC8bl7nxC/pyD1WNk54VqsKLaO9pl6c/2VZwBSykhgevykaG3GDK0j0I5L6fjLQa0NU8/thsQQT6cTtK2Un5YV8rH25ANm7rtB5/lHaVLWla+al6Gcu13K2/a5akgsJuaGfTgWB+Dk3QAm/3mVk3efU9TJitndK9OqovubXWLK4V71ST2z2NjYcPr0aQ4fPsyBAwfo2rUrEydOTFh+4sQJGjRogKOjoW5T586duf7SgIp27dqh0+nw8PDg2TPDEFApJaNGjeLQoUPodDoeP3mMtLTAxsKQ+MMCfQgPfo6NoxPW9um/BBgVHkmwnx9xMREgdFjaOGDn7Jhlzx+8MgEIIe4Q/xTwy6SUxTMlIkVJjUuZl+4JxCcBpxKJVjEz0fF+zSJ0qFKAZf/cZd7BW7wz6zDtvArwxdulKeRolXibvtcNB3+dqeFSk1MJrjwJZuqua+y76oOrrQUT2legS7VCmBnRg0G5kYmJCQ0bNqRhw4ZUrFiR5cuXJyx71Sf1F5eOXl531apV+Pr6cvr0aczMzChatChRUVG4u7sjhMDKzpmIkOeE+D0l7PlzbJ2cyWNrldIuEsRERRPk40dsdBggMM9jj72rEyammXdvITlpuQfw8tNjlkBnIENKHwohlgCtAR8pZYWM2GZmyTZPAL//vuGrMXcFe1OuZQ1JYHlrWNoSCr6V7GpWQH+gdynJHb8w7l0K48olCHHIQwkXGyxM4w/mD44DAj7ayn2Rn+lrvdns/QhbC1OGtSjDx7WLkcc8a/8wlfS7du0aOp2OUqVKAeDt7U2RIkW4eNEwhLh69eoMGTKE58+fY2try8aNG6lYsWKq2wwKCsLV1RUzMzMOHDjAvXuJh/vauThg42hPsF8AkWFBBPk8IvR5HuycnbGwskyyvdiYWIJ8/IiJDAUkZhY22Lk6Y6ZR2Ym0XAL6b+WtGUKII8C3GbD/ZcAc4NcM2Fbu8DCLnn7N7tw8DJU4d3wFAXdSXdUcKGMGxV31+IVGE/Q8mruBAkdrcxytzTFxKkVA48nM/Cea304cRCcEfeuX4NMGJbK+3pDy2kJDQxk4cCCBgYGYmppSsmRJFi5cSKdOnQAoUKAAo0aNokaNGuTPnx8PDw/s7e1T3WaPHj1o06YN1apVw8vLi7JlyyZZR2eiI6+bM3GxDgT5+BMdEczzJw8wNbfG3tUZMwtz4uLiCPbxJyo8GJCYmlth5+KMuaVF0p1mobTUAnr5qV8dhjOCT6WUlTIkACGKAtvScgagZS2g24G3iYiNoLyzNtc3E6iiRG/stm8o0/ZcZ9v5JzhYmfG2hxvbzj8hKlZP17cK8XnjUuSzT/rpTUnZlStXKFeunNZhvFJoaCg2NjbExsbSvn17evXqRfv27TN0HzHRMQT7+hETGQaAqbkVsTERIPWYmFpi6+yMpXWeDN3ny5L7v0ipFlBaLgH99NK/Y4E7QJc3ijAdhBB9gD4AhQsXfsXaGe9ByAN+9v6ZHbd3IJHUzl+bQVUG4eHkkeWxKBmjuIsNc96rQr8GQUzedY11px7S2tOdL5uVoZiztdbhKZlozJgx7N27l8jISJo1a0a7du0yfB9m5mY4FXAnOjKKYF/DdX6diTk2js5Y2WWv36+0nAEUl1Le/s+8YlLK1M+70xpANj0D8IvwY8G5BWy4sQFTYUqPcj3Ia5GXXy7+QlBUEC2KtmBA5QEUsSvy6o1lJHUGkOGiYuOwyOKbbzmNsZwBaEEfp0foRJbVSsroM4ANwH+Lv20AqiazrtELiQ5h6cWlrLyykui4aDqW6kjfSn1xtTLUfu9YuiPLLi1jxeUV7Lm3hw6lOtCvUr+E5ZmuVq1Xr6Okizr4K5kpO5eUTjEBCCHKAuUBeyFEh5cW2WEYDZSjRMZGsubqmoRP+C2LtuSzyp8l+YRva27LwMoD6V62OwvPL2T99fVsubWFHuV60KtCL+wtUr+p9MZ+/DFzt68oSq6R2hlAGQxDNPMCbV6aHwJ8khE7F0KsBhoCzkKIh8B3UsrFGbHttIrVx/LHzT+Ye24uPuE+1ClQh0GVB1HOKfXTWec8zoyqMYoPPD5grvdcll5cyvrr6+lVoRc9yvUgj2nm3eRRFEXJCGm5B1BLSnk0i+JJVUbeA5BSsvf+XmadmcXd4Lt4OnsyuOpg3sqX/JjyV7kWcI3ZZ2fz18O/cMnjQr9K/Whfqj1mugweRtixo+Hrxo0Zu11FeQPqHkD2kZ57AKk1hBkW/8/3hBCz/jtlbMhZ69iTY3Tf3p0vDn6BTuiY0WgGK99Z+doHf4AyjmWY02QOy1ssp6BtQcYdG0e7ze3YeWcnepmBZYr9/Q2ToihpMmvWLMqVK0ePHj2IioqiadOmeHl5sXbt2nRt5+DBg/zzzz+ZFKU2UrsE9KLipzYD7zPBJb9LzDgzg2NPjuFu7c64OuNoU7wNJrqMuwlYxa0Ky1ss5/Cjw8w4M4Nhh4ax5OISBlUZRJ38dXJM1yRFMRZz585l586dFCtWjGPHjhETE4O3t3e6t3Pw4EFsbGyoXbt2xgepkRQTgJRya/zX5SmtYyzuBN1h9tnZ7Lm3BwcLB4a9NYwuZbpgYZI5T+EJIahfsD51C9Rlx50dzDk7h0/3fko1t2qMqD6CMo5lMmW/ipJtvBiu/LIuXaB/fwgPh3feSbq8Z0/D5OcH8U/vJkjjsOdp06axZMkSAHr37s3Vq1e5ffs2bdu25f3332fRokX4+vri5eXFxo0bWbRoEVu2bMHU1JRmzZoxdepUfH196devH/fv3wdgxowZFChQgPnz52NiYsLKlSuZPXs29erVS/PbkV2lNgpoK8kUgXtBStk2UyLKBMsuLePvR3/zaaVP+dDjQ2zM09dYIfzMGXxnzUYfGorzp/2wadw4TZ/kdUJH6+KtaV6kORtubGD+ufn02dOHjW034pzH+XV/HEVRknH69GmWLl3K8ePHkVJSo0YNVq5cyZ9//smBAwdwdnamRo0aTJ06lW3bthEQEMCmTZu4evUqQggCAwMBGDRoEEOGDKFu3brcv3+f5s2bc+XKFfr164eNjQ1fffWVtj9oBkrtEtDULIsik31e+XM+r/w5Tnmc0vW6yGvX8Z0xg9ADBzBxccbEypqHnw0gT6VKuHz5BdbVq6dpO2YmZnQv253q+arTbVs3Rh8Zzbym89CJ1xgf3KRJ+l+jKFkttU/sVlapL3d2fq0HHY8cOUL79u2xtjY8bduhQwcOHz6c4vp2dnZYWlrSu3dvWrVqRevWrQHYu3dvohaSwcHBhITkzHafqV0CSuimIIQwB8piOCO4JqXUtktyOqX3wB/98CF+s2cTtGUrOhsbXIYMwfGD9xHm5gRu2oTfnJ+5/+FHWNerh+sXQ7BM4+iHEnlLMPStoYw7No6Vl1fyYfkP0//DfPNN+l+jKLlAehuzmJqacuLECfbt28eaNWuYM2cO+/fvR6/Xc/ToUfLkyflDuV/5EVQI0Qq4BczCULnzphCiZWYHpoVYPz+ejp/ArZbvEPznLpz+14uSe3bj3LcPjd+xolFTUxw6d6bErj9xHTqUiPPnudO+A4++/Iro/5SJTUnn0p1pXKgx089M54q/6qypKBmlfv36bN68mfDwcMLCwti0aVOq1+lDQ0MJCgrinXfeYcaMGQk3hps1a8acOXMS1nsx39bWNsedCaTlGsRPQCMpZUMpZQOgETmsO1hcaCi+s2Zxs1lznq9eTd727SmxexeuX32FSd68SdbXWVomJAenvn0J2b+fW61a82TsWGJ8fFLdlxCCsbXH4mjhyPDDwwmPCU9fsC1bGiZFURKpUqUKPXv2pHr16tSoUYPevXtTuXLlFNcPCQmhdevWeHp60qBBA6ZPNxzWZs2axalTp/D09MTDw4P58+cD0KZNGzZt2oSXl1eql5aMSVoeBDskpaz/0vcC+OvleVklo4vB6aOieL56Nf7zFxAXGIhtyxa4fP45FsWKJVk3tRpsMT4++M+fz/N16xGmpjh++CFOvf+HiV3KLQiPPznOJ7s/oWPpjnxX67u0B62KwSnZkHoQLPvIkAfBXnJJCLFDCNFTCPERsBU4KYTo8J8aQUZDxsYSuPF3brVoic/ESViWL0/RDRsoOH16sgf/VzFzdSXft99SYsd2bJs2xX/hQm6+3Qz/X35BHxGR7GtquNegV4VebLi+gb339r7pj6QoipJuaUkAlsAzoAGGuj2+GFpCtsFQK8hoSCkJ3rOH2++248no0Zg6O1N42VIKL/6FPBXevNGLeeHCFJg6hWKbfiePVyV8pv7EreYteL52HTI2Nsn6n1X+jApOFfjun+94Gvb0jfevKIqSHmlpCflxVgSS2cKOHcdn2jQiz5/HvHhxCsyaie3bb2fKk7mW5cpReMECwk+exGfadJ5+9x0BS5fiMngQts2aIXSGvGumM2NS/Ul02tqJUUdGsejtRRn6VLKiKEpq0jIKqJgQYpoQ4nchxJYXU1YEl1Gejp/A/Z49ifXxwX3CeIpv+QO7Zs0yvSyD1VtvUeS3VRScOxdhZsqjwUN4+PnniYarFbYrzOgaozn59CRLLy199UZbtzZMiqIobygtDWE2A4sxXPvPwKpmWcemYUPMChTA4b3u6CyytgmzEALbxo2waVAf/0WL8J0xk+crVuL44QcJ67Qt0ZYjj44w5+wcquerjqeLZ8obzEFPISqKoq203AOIlFLOklIekFL+9WLK9MgykE3dOjh93DPLD/4vEyYmOPXti03DhvhMmULktWv/LhOCb2p9g6uVK8MPDScsJkyzOBVFyT3SkgBmCiG+E0LUEkJUeTFlemQ5kBAC9x8moMtrz6Mvv0QfGZmwzM7cjon1JvI47DE/HP8h5Y00bJh8oS1FUdJs2bJlDBgw4I23s379esqVK0ejRo0A6N69O56engnPFKSVt7c3O3bseON40istCaAihg5gEzE8FPYTOahOUFYzdXQk/8SJRN+8hc/kyYmWVXGrQl/Pvmy5tYXtt7drFKGiKGm1ePFi5s6dy4EDB3j69Cn//PMP58+fZ8iQIenajlYJIC33ANoDxY2t/k92ZlOnDo69ehGwZAnWdeti27hxwrI+nn04+vgo44+Np5JLJQraFtQwUkV5Pek9Sc2I5xp//fVXpk6dihACT09PVqxYwdatWxk/fjzR0dE4OTmxatUq3NzcEr3u2bNn9OvXj9u3bwMwb968JDX/V69ezQ8//ICUklatWjFp0iS+//57jhw5wp07d2jbti27du3Cx8cHLy8vZs+ezdmzZ5k/fz6mpqZ4eHiwZs0awsLCGDhwIBcuXCA2NpYxY8bQsmVLvv32WyIiIjhy5AgjR46ka9eub/6GpIWUMtUJWAu4vmq9rJiqVq0qtdKggWHKKPqoKHm7fQd5rUZNGf30WaJlD0Meypqrasoe23vImLiYzA1EUTLA5cuXE33/4tc0rdObunjxoixdurT09fWVUkrp7+8vpZQyICBA6vV6KaWUixYtkl988YWUUsqlS5fKzz77TEopZZcuXeT06dOllFLGxsbKwMDARNt+9OiRLFSokPTx8ZExMTGyUaNGctOmTfE/ZwN58uRJKaWUd+7ckeXLl094nbu7u4yMjJRSSvn8+XMppZQjR46UK1asSJhXqlQpGRoamiieN/Xf/wsppQROyWSOqWk5A3ADrgohTgJR/+YN+W7mpKTcQZibk3/qVO507MjjEcMpvHhxwvMBBWwK8G2tbxl2aBgLzi/gM6/PNI5WUdInqyuV7N+/n06dOuHsbOiz4ejoCMDDhw/p2rUrT548ITo6mmLJPOm/f/9+fv31VwBMTEywt7dPtPzkyZM0bNgQFxcXAHr06MGhQ4do165dqjF5enrSo0cP2rVrl7Du7t272bJlC1OnGq6iR0ZGJjSe0UJa7gF8h+Ey0A/ANOAEUDIzg8otLIoXI9/oUYQfPUbA0sTPALQs1pK2Jdqy8PxCTj87/e+CLl0Mk6IoCaSUyT7XM3DgQAYMGMCFCxdYsGABkS8NvEjPtl/H9u3b+eyzzzh9+jRVq1YlNjYWKSUbN27E29sbb29v7t+/r2kNpVcmAGkY8hkEtAKWAU2A+ZkbVvZz8GDmfKqx79gR2+bN8Zk+g4gLFxMtG1VjFAVsCjDi8AiCooIMM/v3N0yKoiRo0qQJ69atw9/fH4CAgAAAgoKCKFCgAADLlyff3bZJkybMmzcPgLi4OIKDgxMtr1GjBn/99Rd+fn7ExcWxevVqGjRokGo8er2eBw8e0KhRIyZPnkxgYCChoaE0b96c2bNnJySVs2fPAtqVmk4xAQghSgshvhVCXMHQB+ABhuqhjaSUs7MswhxOCIH792MxdXHh8VdfoQ/79xkAazNrJtWbhF+4H+OOjTP80oSHGyZFURKUL1+e0aNH06BBAypVqsQXX3wBwJgxY+jcuTP16tVLuDz0XzNnzuTAgQNUrFiRqlWrcunSpUTL3d3d+fHHH2nUqBGVKlWiSpUqvPtu6lfA4+LieP/996lYsSKVK1dmyJAh5M2bl2+++YaYmBg8PT2pUKEC38Q3eGrUqBGXL1/Gy8uLtWvXZsA7kjYploMWQuiBw8D/pJQ34+fdllIWz7Lo/iOjy0FnJ+EnT3Lvw4+w79Ce/BMmJFr2y4VfmHlmJt/X/p72n8w0zFTloJVsRJWDzj4yqhx0R+ApcEAIsUgI0QTI3OI5uZjVW2/h1K8vQRt/J3jnzkTLPi7/MdXzVefHEz9yzy5OowgVRclpUkwAUspNUsquGHoBHwSGAG5CiHlCiGZZFF+u4tK/P5aVPHny7XfEPHqUMN9EZ8KEuhMwNzFneP0QYnSvd1NKyTwx+hitQ1CUdEvLTeAwKeUqKWVroCDgDYzI7MByI2FmRoGpU0Gv59Gw4ci4fz/t57POx9haY7nkEsc8r+SbzCjaOPTwEO9ufpf7wdoN51OU15GWYaAJpJQBUsoFUsrGr15beR3mhQqR77tviTh9Gr8FCxIta1KkCW1vWrCsQoQ62GQTMXExTDk5BRNhgru1u9bhKEq6pCsBKFnDvm1b7Nq0we/nuYSfOZto2eASPTHTmTHl5BSNolNeturKKu4G32XYW8MwMzHTOhxFSReVALKpfN9+g5m7O4+HDiXupfHBLj0H0LfaAA4+PMjfj/7WMELFL8KP+efnU79gfeoVrKd1OIqSbioBZFMmtrbknzKZmKdPeTpm7L9PI/r58b5LCwrbFmbyycnq5qOGZp6ZSVRcFMPeGqZ1KApw9+5dKlSokGjemDFjEsouKEmpBJCNWVWujMuAzwjevp3gLfFdODt1wrzrewx9ayi3g26z5uoabYPMpS76XWTzzc18UO4DitgV0TocRQNxcXGpfm8M0lIMTtGQU58+hP79N0+/H0eeypUxj5/foGAD6uSvwzzvebQq3gpHS0dN48xNpJT8eOJHnCyd6OPZR+twsp+dI+DphYzdZr6K0HLia7+8YcOG1KhRgwMHDhAYGMjixYupV68ecXFxDB8+nF27diGE4JNPPmHgwIHs27ePr776itjYWN566y3mzZuHhYUFRYsWpVevXuzevZsBAwYwYsSIRN9369YtA3/ozKfOALI5YWJCgcmTwcSER0OHJlwKEkIwrPowImIjmHVmlsZR5i7bbm/jvO95BlUZhI25jdbhKGkUGxvLiRMnmDFjBmPHjgVg4cKF3Llzh7Nnz3L+/Hl69OhBZGQkPXv2ZO3atQl1+1/UCgKwtLTkyJEjCQf7/35vTNQZgBEwy58f9++/59HgwfhaW+NqazjoFLcvTvdy3Vl5eSVdynTBw8lD40hzvvCYcGacnkEFpwq8W1JVRE/WG3xSfxPJVQN9eX6HDh0AqFq1Knfv3gVg79699OvXD1NTw6HQ0dGRc+fOUaxYMUqXLg3ARx99xM8//8zgwYMBkjRrybLmLZlAnQEYCbsWzbHv1BH/sDDCov5tztavUj8cLB2YdGLSa5etVdJu0YVF+ET4MKLGCHTC8OcT88wHv4WL0EdFveLVSmZycnLi+fPnieYFBAQkFIGzsLAADDX/Y2NjgeTLSL/q78ja2jrV742JSgBplB16secbNQozR0ee6gQy2pAE7Mzt+Lzy55zxOcOfd//UNsAc7kHwA5ZfWk6b4m2o5FIpYb7vtJ/wmz2bWB8fDaNTbGxscHd3Z9++fYDh4P/nn39St27dFF/TrFkz5s+fn5AQAgICKFu2LHfv3uXmzZsArFix4pXln42VSgBGRGdlRb5JE4n2DyBgxYqE+e1KtqOcYzl+OvUT4TGqVHRmmXpqKqY6UwZXHZwwL8Lbm6A/tuD48ceYFyqkXXAKYOgLPH78eLy8vGjcuDHfffcdJUqUSHH93r17U7hwYTw9PalUqRK//fYblpaWLF26lM6dO1OxYkV0Oh39+vXLwp8iCyXXJzKrJqAFcA24CYx41fo5qSfwa7t/X97/6CN5tUpVGePjkzD79NPTssKyCnL2mdkaBpdz/fPoH1lhWQW56PyihHn6uDh5u1Nneb1efRkXGqphdNpLrg+too309ATW7AxACGEC/Ay0BDyA7kIIdRfzVT74ALebN9FHR+MzfUbC7CpuVWhZrCXLLi3jUeijlF+vpFuMPoZJJyZR0KYgH3h8kDA/aNNmIi9cwPWrL9EZ8XVgJffS8hJQdeCmlPK2lDIaWAOoYRVpYG5qiuOHHxD0++9EXPh3vPUXVb9AJ3T8dOonDaPLedZdW8etoFsMfWsoFiaGG4lxoaH4TJ9OnkqVsGvdWuMIFeX1aJkACmBoM/nCw/h5iQgh+gghTgkhTvn6+mZZcNmd86efYuLszLPxE5B6PWAoGf2/Cv9jz709nHhyQuMIc4aAyAB+9v6ZWu61aFSoUcJ8v3nziPPzw+3r0QidupWmGCctf3OTG7SbZPyVlHKhlLKalLKai4tLFoRlHExsbHD94gsizp0jeOvWhPkflf+IAjYFmHhyIrH6WA0jzBl+Pvsz4THhDK8+PGG4YNSdOwT8ugL7Dh3IU7GixhEqyuvTMgE8BF4eNlEQeKxRLEbJvt27WHp64jP1p4Rm8pamlnxZ7UtuPL/BhusbNI7QuF0LuMaGGxvoVrYbJfL+O5LEZ+IkdObmuA4ZrF1wipIBtEwAJ4FSQohiQghzoBuwRcN4jMOXXxomQOh05Bs1klhfX/wWLExYpWnhplTPV5053nMIigrSKlKjJuPr/diZ2/FppU8T5oceOkToX3/h3L8/puqMVDFymiUAKWUsMADYBVwB1kkpL2kVj9Fo08Ywxcvj5YX9u20JWLqU6PuGLmFCCIZXH05IdAhzzs7RKlKjtuveLk4/O83AygOxt7AHQEZH8+zHiZgXKYLjB+9rHKHyX0+fPqVbt26UKFECDw8P3nnnHa5fv55h2//hhx8SfV+7du0M27ZWNL17JaXcIaUsLaUsIaWcoGUsRuPaNcP0EpcvvgQzM55Nmpwwr7RDabqU7sK66+u4/jzj/ghyg4jYCH469RNlHcvSsVTHhPkBq34j+s4dXEeOQJibp7IFJatJKWnfvj0NGzbk1q1bXL58mR9++IFnz56l6fVpKe383wTwzz//vH7A2YQqBmds+vY1fD14MGGWmZsrzv364TttGqF//41NnToADKg8gJ13dzL5xGQWNVuUYrEsJbFlF5fxNOwpP9b9EROdCQCxfn74/fwz1vXrYat1TZBsbtKJSVwNuJqh2yzrWJbh1YenuPzAgQOYmZklemLXy8sLMCSHYcOGsXPnToQQfP3113Tt2pWDBw8yduxY3N3d8fb2Zu7cuYm+v3z5csK2RowYQUREBF5eXpQvX55Vq1ZhY2NDaGgoBw8e5LvvvsPNzQ1vb286dOhAxYoVmTlzJhEREWzevJkSJUrg6+tLv379uB9/pj5jxgzqxP+takWNX8shHHt+hFnhwjz78UdkjKFLmL2FPZ95fcbxp8fZd3+fxhEahyehT1hycQnNizanWr5qCfN9ZsxAHxmJ24iRGkanpOTixYtUrVo12WW///473t7enDt3jr179zJ06FCePHkCwIkTJ5gwYULCwf6/378wceJE8uTJg7e3N6tWrUqyj3PnzjFz5kwuXLjAihUruH79OidOnKB3797Mnj0bgEGDBjFkyBBOnjzJxo0b6d27d0a+Ba9FnQHkEDpzc9xGDOdh/894vnoNjh8anljtXLoz66+vZ+qpqdQtUBdLU0uNI83efjpteIjuy6pfJsyLuHiJoI2/49izJxbFi2kVmtFI7ZO6Fo4cOUL37t0xMTHBzc2NBg0acPLkSezs7KhevTrFiv37f/rf79Pqrbfewt3dHYASJUrQrFkzACpWrMiBAwcAQ+nplxNLcHAwISEh2NravsmP90bUGUAOYtOoEda1a+M7Zw6xAQEAmOpMGfHWCB6FPmL5peUaR5i9nXx6kl13d9GrQi/cbQx/zFJKnv3wAyaOjjj3//QVW1C0Ur58eU6fPp3sMplKeeeMKu38otQ0gE6nS/hep9MlVBrV6/UcPXoUb29vvL29efTokaYHf1AJIM0OHkx02T1bEkLgNmok+rAwfGf+2yWsunt13i7yNosvLuZp2FMNI8y+4vRxTDoxCXdrd3pW6JkwP3jbdiLOnMF1yGBMNP5jVVLWuHFjoqKiWLRoUcK8kydP8tdff1G/fn3Wrl1LXFwcvr6+HDp0iOrVq6d7H2ZmZsTEX159Hc2aNWPOnH9H5Xl7e7/2tjKKSgDG5uuvDVMKLEqWxKHHewSuW0fklSsJ87+s9iV6qWf66elZEaXR2XhjI9eeX+PLal+SxzQPAPrwcHymTsWyfHns47tJKdmTEIJNmzaxZ88eSpQoQfny5RkzZgz58+enffv2CeWeGzduzOTJk8mXL1+699GnTx88PT3p0aPHa8U4a9YsTp06haenJx4eHsyfP/+1tpORRGqnR9lNtWrV5KlTp7QOI9uLCw7mVvMWWJQoQeEVvyaM/plzdg4Lzi/g15a/Utm1ssZRZh9BUUG03tSaknlLsqT5koT3y2fGDPznL6DIb79hVUW9X6m5cuUK5cqV0zoMheT/L4QQp6WU1f67rjoDMDbe3oYpFSZ2drgMHkz4qVOE/Plvl7BeFXrhZuXGj8d/JE6fdJxzbjXv3DyCo4MZUX1EwsE/+uFDApYsxa5NG3XwV3IslQCMzeDBhukV8nbqiEW5cjybPAV9RAQAVmZWfFH1C64EXGHjjY2ZG6eRuBZwjTVX19CpVCfKOJZJmO8zaTKYmOD65RcaRqcomUslgBxKmJiQb/QoYp88wf+XxQnzWxZrSY18NZh6aip3g+5qF2A2EBkbyYjDI8hrkZcBlQckzA87doyQPXtw7tsHs9e4VqwoxkIlgBzMqlo17N5pif8vvxDzyNAlTAjB+LrjMTcxZ/jh4cTEvf6oBmM37fQ0bgbeZHzd8ThYOgAgY2N5NuEHzAoWxPHjjzWOUFEyl0oAOZzr0KEgBM+mTk2Yl886H2NrjeWy/2Vme8/WMDrt/PXgL1ZfXc0HHh9Qt0DdhPnP16wl6sYNXIcPQ/fS2G5FyYlUAsjhzNzdcfqkNyE7/yTsxL9dwpoUaULn0p1ZenEpRx8f1TDCrOcb7ss3f39DGYcyDK4yOGF+7PPn+M6ejVWtmtg2bapdgIqSRVQCMDY//GCY0sGpVy9M87vz7IcfkS9VORz61lCK2Rdj9JHRPI98ntGRZkt6qWf0kdFExEYwuf5kzE3+rerpN3s2+tBQ3EaOVIXzjFBml4NOScOGDTHW4ekqARib2rUNUzro8uTBbdgwoq5eJXD9+oT5eUzzMLn+ZAKjAvn2n29TfWQ+p1hxeQVHnxxlWPVhFM9bPGF+5LVrPF+zFodu3bAsXVrDCJXXkRXloHMiVQzO2LyoQZ7OJGDbvDlWb72F74yZ2LVsiYm9oclJWceyDKk6hMknJ7Pu2jq6lu2a0RFnG5f9LzPjzAyaFG5Cp1KdEuZLKXk24QfD8xMDB6SyBSUtnv7wA1FXMrYctEW5suQbNSrF5ZldDhrAxsaGvn37cuDAARwcHFizZg0v+pSvX7+e/v37ExgYyOLFi6lXrx53797lgw8+ICy+XeucOXOoXbs2T548oWvXrgQHBxMbG8u8efOoV68eu3fv5rvvviMqKooSJUqwdOlSbGxsMvR9/C91BmBsRo0yTOkkhMDt69HEBQfjO+fnRMt6lOtBnQJ1mHJqCrcCb2VUpNlKeEw4ww8Nx9HSkTG1xiS6xBOyew/hJ07gMuhzTPLm1S5I5bVldjlogLCwMKpUqcKZM2do0KABY8eOTVgWGxvLiRMnmDFjRsJ8V1dX9uzZw5kzZ1i7di2ff/45AL/99hvNmzdPiMnLyws/Pz/Gjx/P3r17OXPmDNWqVWPatGkZ+h4lR50B5CKWZcqQt2sXnv/2Gw5dOmNRqhQAOqFjfJ3xdNzSkWGHhvFbq9+wMMlZI2Amn5zMveB7/NLsF/Ja5k2Yr4+MxGfSJCzKlCFvly7aBZiDpPZJXQsZVQ5ap9PRtavhDPn999+nw0v1oV78u2rVqty9exeAmJgYBgwYgLe3NyYmJgn3I9566y169epFTEwM7dq1w8vLi7/++ovLly8nNIiJjo6mVq1aGf5eJPmZMn0PSrbi8vnn6GxseDz664TGMQDOeZwZV2cc159fZ8bpGdoFmAl2393Nxhsb+V/F/1HdPXEVSN/p04l5/Bi3UaMQJiYaRai8KS3KQb98Fvmi/LOJiUlC+efp06fj5ubGuXPnOHXqFNHR0QDUr1+fQ4cOUaBAAT744AN+/fVXpJS8/fbbCaWiL1++zOLFi5PuNIOpBJBGDRsaJmNn6uCA+9ixRJ4/n+RSUP2C9Xm/3PusvLKSQw8PaRRhxnoa9pQxR8dQwakC/b36J1oWeugQAct/xeH997Gukf7ywEr2kRXloPV6PRs2bAAMl3Hq1q2b6vpBQUG4u7uj0+lYsWJFwo3le/fu4erqyieffML//vc/zpw5Q82aNfn777+5efMmAOHh4VkygkklgFzIrkVz7Dt1xH/hQsKOn0i0bHDVwZR2KM03f3+DX4SfRhFmjDh9HCMOjzDU+q8/CTOdWcKyWD8/Ho8chUWpUrgO/UrDKJWMkBXloK2trbl06RJVq1Zl//79fPvtt6mu379/f5YvX07NmjW5fv16wtnFwYMH8fLyonLlymzcuJFBgwbh4uLCsmXL6N69O56entSsWZOrVzP2RnpyVDnoNHrx6V/zpjAvKoHGj3B4XfqwMO507IQ+IoLif2xOdPPzVuAtum7rSrV81ZjbZC46YZyfExaeX8jss7MZX2c875Z8N2G+1Ot50Lcf4SdOUHT9OjXsMwPkhnLQL5rAZ3eqHHRO5uX1xgd/AJ21NfmnTiU2IIAn3yR+BqBE3hIMrTaUvx/9zaorSRtgG4NzvueY6z2XlkVb0rZE20TLnq9cSdjhw7gOH6YO/kquphKAsdm71zBlgDwVyuM6eDAhe/YQGH9t84UuZbrQsFBDpp+eztWAzD8VzUih0aEMPzQcNys3vq71daKbdZFXr+IzZSo2jRrh0L27hlEqxsYYPv2nl0oAxmb8eMOUQRw/7ol17Vo8++FHom7fTpgvhOD72t+T1yIvww4NIyI2IsP2mdkmHJ/Ak7AnTKw/ETtzu4T5+ogIHn35FSZ58+L+wwRV7kHJ9VQCyOWETof7xInoLC159NVX6OOHqgE4WDrwQ70fuBt0lyknp2gYZdptvbWVbbe30c+zX5K2l88mTSL61i3yT5qIqYODRhEqSvahEoCCmasr7hMmEHX5Cr7TZyRaVtO9Jj0r9GT99fXsu7dPmwDT6EHIAyYcn0Bl18p84vlJomUhe/cSuGYtjv/rhXU6y2goSk6lEoACgG3jRji89x4BS5cSeuTvRMsGeg3Ew8mD745+x9OwpxpFmLoYfQwjDo1Ah46J9SZiqvv3IfeYZ894MvprLD08cB00SMMoFSV7UQlASeA6bCgWpUryeMQIYv39E+abmZgxqd4kouOiGX1kdLZsKD//3HzO+53nm1rfkN8mf8J8GRfH42HD0UdHk/+nqQhz81S2ohizlMpBHzx4kNatW2fIPjZv3pxsnaDkln377bfszaABG5lFJQBjs2CBYcoEOktL8k/9CX1wME9GjU40NLSofVFGVh/JiacnWHppaabs/3WdfHqSRecX8W6Jd2lZrGWiZf5LlhB+/Dj5vh6NRQo1XhTj96bloNMqPQng+++/p2k2byykisEZmzJlMnXzlmVK4zpsGM/Gj+f5qt9wfL9HwrJ2Jdvx9+O/+fnsz9TIV4OKLhUzNZa0CIoKYuThkRSyLcTIGiMTLYu4cAHfmbOwbdEC+5cKdymZ6/C66/g9yNghk86FbKjXJeVnNlIrB33w4EFCQ0Pp1KlTQtXQlStXIoTg9OnTfPHFF4SGhuLs7MyyZctwd3dn0aJFLFy4kOjoaEqWLMmKFSvw9vZmy5Yt/PXXX4wfP56NGzdSokQJAP75558ky8aNG0fr1q3p1KkTRYsW5b333uPAgQPExMSwcOFCRo4cyc2bNxk6dGhC3FOmTGHdunVERUXRvn37RBVHM4M6AzA2W7capkzk0OM9bBo0wGfyZCKv/VuPRAjBNzW/wcXKheGHhxMare24aCklY4+OxT/Cn0n1J2Ft9m8hr7jQMB599RWmLi64jx2jhnzmcKmVgwY4e/YsM2bM4PLly9y+fZu///6bmJgYBg4cyIYNGzh9+jS9evVi9OjRgKG658mTJzl37hzlypVj8eLF1K5dm7Zt2zJlyhS8vb0TDv5AqsteKFSoEEePHqVevXr07NmTDRs2cOzYsYSSErt37+bGjRucOHECb29vTp8+zaFDmVuTS50BGJuffjJ8bdMm03YhhMD9xx+4/e67PP7qS4quX4/O0hIAewt7JtabyMe7PqbdH+3o79WftiXaJrrpmhWOPznOjNMzuOh/kUFVBlHBuUKi5c8mTCDmwUOK/Lo8ofmNkjVS+6SulerVq1OwYEHAcGZw9+5d8ubNy8WLF3n77bcBQxcwd3d3wJBQvv76awIDAwkNDaV58+ZvHEPbtoYn0itWrEhoaCi2trbY2tpiaWlJYGAgu3fvZvfu3VSubBi+HBoayo0bN6hfv/4b7zslKgEoyTJ1dCT/jxN50Ls3PpOnkO/bbxKWVXGrwpLmS5h2ahrf/fMdSy8u5fMqn9O0cNNM/6R9ye8SM8/M5OiTo+Szzsf3tb+nXcl2idYJ2r6doE2bcPq0H1bVkpQ/UXKg8uXLJ1TqTM6Lcs3wb8lmKSXly5fn6NGjSdbv2bMnmzdvplKlSixbtoyDGVAE7EUMOp0uUTw6nS4hnpEjR9K3b9833ldaqUtAaXTwYDYoBJfFbOrWwbFnT57/9hsh+w8kWlbVrSor31nJjEYz0AkdXxz8gve2v8exJ8cyJZY7QXf48uCXdNvejSsBVxhabSjb2m+jfan2iZJO9MNHPB0zljyVKuHSv38qW1RyktTKQaekTJky+Pr6JiSAmJgYLl26BEBISAju7u7ExMSwatW/9bBsbW0JCQlJdnupLUuL5s2bs2TJkoSSE48ePcLHx+e1t5cWKgEoqXL5YggW5crxZNQoYp4l/mUUQtCkcBN+b/s74+qMwy/Sj092f0Kf3X245HcpQ/b/NOwpY/4ZQ/s/2nP40WH6VerHzg47+bD8h0m6lsnYWB4PGwZ6PfmnTkGYmaWwVSWnSa0cdErMzc3ZsGEDw4cPp1KlSnh5efFPfM/tcePGUaNGDd5++23Kli2b8Jpu3boxZcoUKleuzK1bidunprYsLZo1a8Z7771HrVq1qFixIp06dXqjhJIWqhy0sdGgLnXU7dvc6dARqyqVKfTLLwhd8p8bouKiWHdtHYvOL+J51HOaFWnGgMoDKGaf/uGXQVFBLL6wmN+u/kacjKNrma58UvETnPI4pfga3zk/4zdnDvmnTMY+E++RKEnlhnLQxiI95aDVPQBjs2JFlu/Sonhx3EaN5Om33xGwdBlO/+uV/HomFnzg8QHtS7bn18u/svzScvbd30e7ku3oV6kf+axf3YQjPCacVVdWsfTiUkJjQmlTog39vfpTwKZA6q87cwa/uXOxa9tGHfwVJY00OQMQQnQGxgDlgOpSyjR9rFdnANqRUvLo80GEHDxI0dWryVOh/Ctf4x/hz6ILi1h7bS0mwoT3yr7H/yr+D3uLpKNyYuJi2HhjI/PPzcc/0p+GhRoysPJASju8ekRJXHAwd9q1B52OYps3YWJj81o/o/L61BlA9mEMDWEuAh2AnNF4NiutXWuYspgQAvdx32Pq6MjjL79EHxb2ytc45XFiRPURbGu/jeZFm7Ps0jJabmzJovOLCI8JB0Av9Wy/vZ22m9sy4fgEitgV4deWvzK78ew0HfyllDwdM4aYZ88oMHWKOvgrSjpocglISnkFUA/nvI558wxfu3bN8l2b5M1L/smTud+zJ09//JH8aexLUMCmABPqTqBn+Z7MOjuLWWdnserKKrqW7cq+e/u49vwapR1K83OTn6lXoF66fi+CNv9B8I6duAweRJ4M6JSmKLmJugegpIt1jeo49emD/4IFmDo44tS3T5o/dZdyKMXsxrPx9vFm+unpzPWeS0GbgkysN5GWxVqmq/ewjI0l6I8tPJswAatq1XD65JNXv0hRlEQyLQEIIfYCyd31Gy2l/CMd2+kD9AEoXLhwBkWnvAmXAZ8R8/gx/osWEbh+PU79+uLQvTs6C4tXvxjwcvViWYtl3A+5T37r/JiZpH24ppSS0H378Jk+g+hbt7CsWNEw5NPE5HV/HEXJtTLtHoCUsqmUskIyU5oP/vHbWSilrCalrObi4pJZ4SrpIMzMKDBlMkU3bMCyfHl8Jk7iVouWBG78HRkbm7ZtCEERuyLpOviHHT/B3W7deDhgIOj1FJg1k6Lr1mKW79Wji5ScT+ty0ClZtmwZAwYMyJD9ZzT1IJjy2vJUKE/hxb9QeNlSTJ2deTJ6NLffbUfI3r1k5OiyiEuXuN/7E+5/9BGxz3xwHz+O4lu3YNesmbqPpADZoxy0MdLkHoAQoj0wG3ABtgshvKWUb15tKTdIpd6JVqxr1qTourWE7NmD74yZPBwwEMtKnrh+8SXWNaq/9naj797Fd9YsgnfsxMTeHtdhw3B4r3tCYTolezqwbCE+925n6DZdixSnUc8+Ke9T43LQYKgfZGlpyaVLl3j27BnTpk1LOPN4/PgxLVq04NatW7Rv357JkycD8Omnn3Ly5EkiIiLo1KlTQvnnESNGsGXLFkxNTWnWrBlTp07F19eXfv36cf/+fQBmzJhBnTp13uh91WoU0CZgkxb7NnrOzlpHkCwhBHbNmmHbuDFBf/yB7+w53P/oI6zr1sVlyGDylH/1cwMvxDzzwW/uXAI3bECYm+P0aT+cevXCxNY2E38CxZilpRz0pUuXyJ8/P3Xq1OHvv/+mRo0aDBw4kD/++AMXFxfWrl3L6NGjWbJkCR06dOCT+IEFX3/9NYsXL2bgwIG0bds2ocZ/cu7evctff/3FrVu3aNSoETdv3gTA29ubs2fPYmFhQZkyZRg4cCCFChViwoQJODo6EhcXR5MmTTh//jwFCxZk06ZNXL16FSEEgYGBAAwaNIghQ4ZQt25d7t+/T/Pmzbly5cobvW9qFJCxWbbM8LVnTy2jSJEwNSVvx47YtW7N899W4z9/Pnc7dsLunZa4fP455kWLpvjauKAg/H/5hYAVK5FxcTh064Zzv76Yqns/RiW1T+payapy0F26dEGn01GqVCmKFy/O1atXAWjSpAn28WXJPTw8uHfvHoUKFWLdunUsXLiQ2NhYnjx5wuXLl/Hw8MDS0pLevXvTqlWrhLOIvXv3Jrr8FBwcTEhICLZv8MFIJQBjk80TwAs6CwucPu5J3k4d8V+yhIBlywnetZu8nTrh3L8/Zm6uCevqIyIIWLkS/0W/oA8Jwa51a1w+H4h5oUIa/gSKMcku5aD/e0/qxffJ7f/OnTtMnTqVkydP4uDgQM+ePYmMjMTU1JQTJ06wb98+1qxZw5w5c9i/fz96vZ6jR4+SJ0+eNMWSFuomcBo1bPhvHTYl7UxsbXEdNIiSe3bj0K0bgb//zq3mzfH56Sdi/f15vmYNt5o1x/enaVhVrkyxTb9TYMpkdfBX0iU7lIMGWL9+PXq9nlu3bnH79m3KpNLCNTg4GGtra+zt7Xn27Bk7d+4EDI1ggoKCeOedd5gxYwbe3t6AoVronDlzEl7/Yv6bUAlAyRKmzs7k++ZrSuzYjm2zt/H/ZTE36tbj6ZixmBUqRJGVKyi0YD6WL5XeVZS0yg7loMGQVBo0aEDLli2ZP38+lqkMWKhUqRKVK1emfPny9OrVK+GGbkhICK1bt8bT05MGDRowffp0AGbNmsWpU6fw9PTEw8OD+fPnv9Z79TJVDjqNNKjCnLxsE8ibibx2jaDff8eqZk1sGjZUwzmNnCoGZ7hslNoN4qyiykEr2Z5lmTJYjhypdRiKkqupBGBsduzQOgJFUZKx7MUADSOiEoCxsbLSOgJFSZaUUl3K01h6L+mrm8DGZu5cw6Qo2YilpSX+/v4ZWgJESR8pJf7+/qneeP4vdQZgbNatM3zt31/bOBTlJQULFuThw4f4+vpqHUquZmlpmfDAW1qoBKAoyhszMzOjWLFiWoehpJO6BKQoipJLqQSgKIqSS6kEoCiKkksZ1ZPAQghf4N5rvtwZ8MvAcDKLscQJxhOrijPjGUusKk6DIlLKJGV1jSoBvAkhxKnkHoXObowlTjCeWFWcGc9YYlVxpk5dAlIURcmlVAJQFEXJpXJTAliodQBpZCxxgvHEquLMeMYSq4ozFbnmHoCiKIqSWG46A1AURVFeohKAoihKLpUrEoAQooUQ4poQ4qYQYoTW8aRECHFXCHFBCOEthNCm9VkyhBBLhBA+QoiLL81zFELsEULciP/qoGWML6QQ6xghxKP499VbCPGOljHGx1RICHFACHFFCHFJCDEofn62el9TiTNbvadCCEshxAkhxLn4OMfGz89u72dKcWryfub4ewBCCBPgOvA28BA4CXSXUl7WNLBkCCHuAtWklNnqwRUhRH0gFPhVSlkhft5kIEBKOTE+qTpIKYdrGWd8XMnFOgYIlVJO1TK2lwkh3AF3KeUZIYQtcBpoB/QkG72vqcTZhWz0ngpDIwJrKWWoEMIMOAIMAjqQvd7PlOJsgQbvZ244A6gO3JRS3pZSRgNrgHc1jsmoSCkPAQH/mf0usDz+38sxHBQ0l0Ks2Y6U8omU8kz8v0OAK0ABstn7mkqc2Yo0CI3/1ix+kmS/9zOlODWRGxJAAeDBS98/JBv+AseTwG4hxGkhRB+tg3kFNynlEzAcJABXjeN5lQFCiPPxl4iyxeWqF4QQRYHKwHGy8fv6nzghm72nQggTIYQ34APskVJmy/czhThBg/czNySA5HrUZdfrXnWklFWAlsBn8ZczlDc3DygBeAFPgJ80jeYlQggbYCMwWEoZrHU8KUkmzmz3nkop46SUXkBBoLoQooLGISUrhTg1eT9zQwJ4CBR66fuCwGONYkmVlPJx/FcfYBOGy1fZ1bP468MvrhP7aBxPiqSUz+L/6PTAIrLJ+xp/DXgjsEpK+Xv87Gz3viYXZ3Z9TwGklIHAQQzX1bPd+/nCy3Fq9X7mhgRwEiglhCgmhDAHugFbNI4pCSGEdfxNNoQQ1kAz4GLqr9LUFuCj+H9/BPyhYSypenEAiNeebPC+xt8MXAxckVJOe2lRtnpfU4ozu72nQggXIUTe+H/nAZoCV8l+72eycWr1fub4UUAA8UOqZgAmwBIp5QRtI0pKCFEcw6d+MLTq/C27xCmEWA00xFCy9hnwHbAZWAcUBu4DnaWUmt98TSHWhhhOrSVwF+j74rqwVoQQdYHDwAVAHz97FIbr69nmfU0lzu5ko/dUCOGJ4SavCYYPtuuklN8LIZzIXu9nSnGuQIP3M1ckAEVRFCWp3HAJSFEURUmGSgCKoii5lEoAiqIouZRKAIqiKLmUSgCKoii5lEoAiqIouZRKAEquIYRweqnc7tOXyu+GCiHmZsL+lgkh7ggh+r3m6w/Ex1Yto2NTFDA8cKQouYKU0h/DwzZZWSJ6qJRyw+u8UErZSAhxMIPjUZQE6gxAyfWEEA2FENvi/z1GCLFcCLFbGBr0dBBCTBaGRj1/xtfFQQhRVQjxV3zl1l3/eZQ/pf0sE0LMEkL8I4S4LYToFD/fXQhxKP5s5KIQol7m/sSKYqASgKIkVQJohaGW/ErggJSyIhABtIpPArOBTlLKqsASIK1lO9yBukBrYGL8vPeAXfEVIisB3hnzYyhK6tQlIEVJaqeUMkYIcQFDzZY/4+dfAIoCZYAKwB5DrTRMMJTwTYvN8RUfLwsh3OLnnQSWxCeWzVJK7wz5KRTlFdQZgKIkFQUQf6COkf8WzNJj+NAkgEtSSq/4qaKUsll6th1PxO/nEFAfeASsEEJ8mBE/hKK8ikoAipJ+1wAXIUQtMNTLF0KUf92NCSGKAD5SykUYSi9XyZgwFSV16hKQoqSTlDI6/gbuLCGEPYa/oxnApdfcZENgqBAiBkNDe3UGoGQJVQ5aUTKJEGIZsO11h4HGb+Mg8JWU8lRGxaUoL6hLQIqSeYKAcW/yIBhQHIjJ0KgUJZ46A1AURcml1BmAoihKLqUSgKIoSi6lEoCiKEoupRKAoihKLvV/01neTJrm13YAAAAASUVORK5CYII=\n",
|
|
"text/plain": [
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"#####\n",
|
|
"# Plot signal\n",
|
|
"#####\n",
|
|
"\n",
|
|
"fig, ax = plt.subplots()\n",
|
|
"ax.set_xlabel(\"Time [ns]\")\n",
|
|
"ax.set_ylabel(\"Amplitude\")\n",
|
|
"\n",
|
|
"\n",
|
|
"ax.plot(sig1_time/ns, sig1+2, label='Signal')\n",
|
|
"ax.axvline(sig1_offset/ns, color='r', ls='--', label='offset')\n",
|
|
"\n",
|
|
"\n",
|
|
"ax.axvline(calc_lag_time/ns, color='b', ls=(2, (10, 10)), label='calc offset')\n",
|
|
"ax.plot(corr_time/ns, corr_sig+2, label=\"Uncorr\")\n",
|
|
"ax.plot(sig1_time/ns, sin_delay(f_sine, sig1_time, t_delay=calc_lag_time), label=\"Corr time\")\n",
|
|
"ax.plot(sig1_time/ns, sin_delay(f_sine, sig1_time, phase=-1*calc_phase_lag)-0.2, label=\"Corr phase\")\n",
|
|
"\n",
|
|
"# debugging\n",
|
|
"if True:\n",
|
|
" ax.plot(corr_time[-1]/ns + corr_time/ns, sin_delay(f_sine, corr_time, t_delay=sig1_offset)+2, label=\"Cheat time\")\n",
|
|
" ax.plot(corr_time[-1]/ns + corr_time/ns, sin_delay(f_sine, corr_time, phase=-2*np.pi*sig1_offset*f_sine)+2, label=\"Cheat phase\")\n",
|
|
" \n",
|
|
" \n",
|
|
"ax.legend()\n",
|
|
"fig.show();\n",
|
|
"print((sig1_offset - calc_lag_time)*samplerate)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 8,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stderr",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"<ipython-input-8-b80dbc5f3c93>:20: UserWarning: Matplotlib is currently using module://ipykernel.pylab.backend_inline, which is a non-GUI backend, so cannot show the figure.\n",
|
|
" fig.show();\n"
|
|
]
|
|
},
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEHCAYAAACncpHfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABEHUlEQVR4nO3de1xUZf7A8c/DHQS5DQgC5gVFBRQLU9ESy8pMLctL22Wjtqy1tmxrt7LLVpbttrZru5Wrv61oWytdzVI3LTEva6KJCYIKXvCGyF0R5A7P748ZDAFhgJk5M8zzfr3Oi5k5Z8759jTOd85znvN9hJQSRVEUxf44aB2AoiiKog2VABRFUeyUSgCKoih2SiUARVEUO6USgKIoip1SCUBRFMVOOWl5cCHECaAMqAfqpJSxbW2v0+lk3759LRCZ0q6sLP3fiAht47ACWcX6tojw73xbVBRXAODh72GSmBSlqb179xZJKQOav65pAjCYIKUsMmbDvn37kpKSYu54FGPEx+v/bt2qZRRWIT4xHoCtCVs7vY/E+EQAErYmdDkeRWlOCHGytddVF5CiKIqd0voMQALfCSEksFRKuUzjeBRjvfSS1hFYjZeuV22h2CatE8BYKWWuECIQ2CSEyJRSbm+6gRBiDjAHoE+fPlrEqLRm4kStI7AaE/urtlBsk6ZdQFLKXMPfAmANcG0r2yyTUsZKKWMDAlpcw1C0kpqqXxRS81JJzUvVOgxF6TDNzgCEED0AByllmeHxzcDrWsWjdNC8efq/6iIw8zbOA7p2EVhRtKBlF1AvYI0QojGOz6SUGzWMR1EUxa5olgCklNnAcK2OryjWRA3/VLSg9UXgDskqzro05rrRrMhZzB05l4raCiYvn9ziPQkxCSTEJFBUUcSMlTNarP917K+ZHTWb06WnuX/N/S3WPzPmGaZGTCWrKItH1z/aYv1L17/ExP4TSc1LvdQV0NTCGxcSFxbHztM7mb95fov1iyctJiYohqTsJN7Y/kaL9UunLCVCF8G6rHW8k/xOi/WfTv+UMO8wVmSsYEnKkhbrV81ahc5DR2JqIompiS3Wf3PvN3g4e/DBng9YeWBli/WN3RqLdi5i/eH1P8edl4qDcGCY4fmCbQvYfHzzZe/19/Bn9azVALyQ9ALJOcmXrQ/tGcq/7/w3oO9Gad6PPsh/EMum6geGzVk3h8PFhy9bHxMUw+JJiwG478v7yLmQc9n6MaFjeGviWwDctfIuiiuKL1t/Y78beXn8ywDcuvxWKmsrL1s/ZdAUno17FqDF5w5+/uxV19aRlr+f2KXj8HT9+Z+U+uyZ57MH4O7szoZ7NwD2/dnr7PdeI5tKAIpijc5eqKSqrp6MM6U4Ozrg4+GMj4cLFTX1WoemKG0StjQjWGxsrFR3AluJnTv1f+PitI1DY1JKhr31HkE9XZk7djJbMgvZdriQ0spaHB0E11zlyw2DA5kQEcigXp4YrnkpikUJIfa2VmpHJQBF6YK9J89x15KdvDNzOHddEwpAXX0DqafPsyWrgC2ZhRw8ewGA3t5uxA8O5IaIQOLC/fFwUSfgimVcKQGoT6DSOeoMAIB1abnUO2fi7e0F6BOAk6MDsX39iO3rx+9uGUxeaRVbswrYklXA1/vO8NnuU7g4OjCqv9+ls4O+uh7a/ocodkmdASido4rBUd8gGf3WZs44PcegXl5G3QdQU9fAnhMlbMnUJ4RjhRcB6K/rwYzYUObGh5s5asUeqTMARTGx3dnFFJZV4x/mavR7XJwcGBuuY2y4jpemDOVUcQVbsgr43z9S+HLLSSZEBDIkuKcZo1aUn6lqoIrSSev25+Lh4oivh3On99HH34MH4voy+ngZAzOK+WrfGRNGqChtUwlAUTqhpq6BDRl53DS0Fw4mGNnj5CjwcXfm69RcGhpsp1tWsW0qAShKJ/xwtIjzFbVMHdbbZPvUebqSd6GKXceL299YUUxAXQNQOmfxYq0j0NS6tFx6ujlx/aAAFvstNsk+fXs44+nqxFf7zhA3QGeSfSpKW9QZgNI5MTH6xQ5V1dbz7YE8bo0KxsXJgZigGGKCYrq8XwchuCUyiA3peVTVqruIFfNTCUDpnKQk/WKHtmQWcLGmnqnD9d0/SdlJJGWbpi2mjwihrLqO7zMLTLI/RWmL6gJSOucNQ/EwO5wZbN3+XHSeLozu7wdwqZCaKWYGGzPAn0AvV9bsO8Pk6OAu709R2qISgKJ0QHl1HZsPFTB7ZBhOjqY7gW5aDvr2mN4k7jzB+YoafDxcTHYMRWlOdQEpSgckHcynuq6BacNNN/qnudtjQqitl/w3/azZjqEooBKAonTIurRcenu7cXUfX7MdI7J3TwYGeqqbwhSzUwlAUYx0vqKG7UcKmTK8Nw4O5ivrLITgjhEh7DlxjtMlFWY7jqJofg1ACOEIpABnpJRTtI5HMdLSpVpHYHEbM/KorZctbv5aOsX0bXF7TG/+/G0Wa9NyeXyCKhCnmIc1nAE8BRzSOgilgyIi9IsdWbc/l77+HkSFXF6sLUIXQYTOtG0R6uvBtX39+PKnHGypYq9iWzRNAEKIUOA24J9axqF0wrp1+sVOFJRVkXysmGnDe7eY1Wtd1jrWZZm+Le4YEcKxwoscyL1g8n0rCmh/BrAY+D3QoHEcSke9845+sRMb0vNokFy6+aupd5LfaXXS9I5ITUwlNTH1stcmRwfh7CjUxWDFbDRLAEKIKUCBlHJvO9vNEUKkCCFSCgsLLRSdolxubVoug4O8GNjLyyz7by0B+Hi4MCEikK/TcqlXFUIVM9DyDGAsME0IcQL4ArhBCPHv5htJKZdJKWOllLEBAQGWjlFRyDlXwd6T51r99W9u00eEUFhWzc5jRRY/ttL9aZYApJQvSClDpZR9gbuB76WU92kVj6JcyX/362/ImjLM8qUZJgwOxMvNia/25Vr82Er3p/U1AEWxeuv25zI8zIer/C0/cbubsyOTo4LZmHGWyhpVIVQxLc3vAwCQUm4FtmochtIRn36qdQQWkV1YTsaZC7x025ArbvPpdPO2xR0jQliRcppNh/LNWoJCsT/qDEDpnLAw/dLNrd9/FiFgShszf4V5hxHmbb62GNXPj2BvN75Wo4EUE1MJQOmcFSv0SzcmpWRtWi4j+/oR5O12xe1WZKxgRYb52sLBQTAtpjfbDhdSXF5ttuMo9kclAKVzlizRL91YZl4ZRwvK2x39syRlCUtSutYWCVsTLisJ3dz0ESHUNagKoYppqQSgKFewLi0XRwfB5KggrUNhcFBPBgd5sUZ1AykmpBKAorRCSsm6/bmMDdfh7+mqdTiA/mLwvlPnOVl8UetQlG5CJQBFaUVaTimnSyqZqsHY/yvR1yFC3ROgmIxKAIrSirWpubg4OnBzpPbdP416+7gzup8/X6WeURVCFZOwivsAFBu0apXWEZhNfYNk/f5cxkcE4O3u3O72q2ZZri3uGNGb51ansz+nlOFhPhY7rtI9qTMApXN0Ov3SDe05UUJBWbXRN13pPHToPCzTFpOignFxclAXgxWTUAlA6ZzERP3SDa1Ly8Xd2ZEbhwQatX1iaiKJqYldOmZr1UBb4+3uzMQhgazfn0tdvaqirnSNSgBK53TTBFBb38CGjDwmDu2Fh4txPaSWTAAAt8eEUFRew46jqkKo0jUqAShKEzuPFVNyscaqRv80F2+4NqEmilG6SiUARWlibWouXm5OjI+w3rknXJ0cuW1YMN8eyOdidZ3W4Sg2TCUARTGoqq3nuwN5TIoMwtXJUetw2nRHTAiVtfV8dzBP61AUG6YSgKIYbDtcSFl1nSYzf3VU7FW+hPi4q5vClC5R9wEonfPNN1pHYHLr0nLx6+FC3AD/Dr3vm3st3xYODoI7RvRmydZjFJZVE+BlHeUqFNuizgCUzvHw0C/dREVNHZsPFTA5Oggnx479s/Bw9sDD2fJtcUdMCA1Sn7gUpTPUGYDSOR98oP87d662cZjIpoP5VNbWM7WNiV+u5IM9+raYO7LzbdFWKegrGdjLi8jePfk69QwPjevX6WMr9kuzMwAhhJsQ4kchRJoQ4oAQ4jWtYlE6YeVK/dJNrEs7S1BPN0b29evwe1ceWMnKA9q0xfQRIaTllHKssFyT4yu2TcsuoGrgBinlcCAGmCSEGK1hPIqdKq2oZdvhAqYMC8bBQWgdTodMHd4bB4GaLlLpFM0SgNRr/NnibFhUiUPF4r49mEdtvbSJ0T/N9erpRtwAHV+l5qoKoUqHaXoRWAjhKIRIBQqATVLK3VrGo9indWm59PHzYFiot9ahdModI0I4VVLBT6fOax2KYmM0TQBSynopZQwQClwrhIhqvo0QYo4QIkUIkVJYWGjxGJXurai8mp3Hipk6PBghbKv7p9Etkb1wc3ZQpSGUDrOKYaBSyvPAVmBSK+uWSSljpZSxAQHWe3u+3dm6Vb/YuO8zC6hvkNwW3fnun60JW9masNV0QXWQl5szE4f0Yv3+XGpVhVClA7QcBRQghPAxPHYHJgKZWsWj2Kddx4rx7+HCkGAvTePoSDXQ1kwfEcK5ilq2H1ZnyYrxtDwDCAa2CCH2A3vQXwNYr2E8SkcsWqRfbJiUkuTsYkb39+9S98+inYtYtLNrbdHVBHD9oAB8PZzVRDFKh2h2I5iUcj8wQqvjK1203pCrn31W2zi64GRxBWdLqxjdwdIPza0/rG+LZ+O0awtnRwcmRwfz5U9nqK6rt/pidop1sIprAIqiheTsYgDG9O9aArAW8RGBVNbWs0+NBlKMpBKAYreSjxUT4OXKgIAeWodiEqP6++HoIPhBzRSmGEklAMUuNfb/j+li/7816enmzPBQbzVVpGI0lQCUznF31y826ljhRQrLqhnTxf5/AHdnd9ydraMtxoXrSDt9ntLKWq1DUWyAqgaqdM6GDVpH0CWm7P/fcK/1tMXYcB1/+/4ou7KLuSUySOtwFCunEoBil3YdKybY242r/K1jToPOlINuzYg+vrg7O/LD0SKVAJR2qS4gpXMWLNAvNkhKyS4T9v8v2LaABdusoy1cnBwY1d9PXQdQjKISgNI5mzfrFxt0OL+c4os1XR7/32jz8c1sPm49bTEuXEd24UVyz1dqHYpi5VQCUOxO8jH9r+PuMv6/uXEDdQBqOKjSLpUAFLuz81gxob7uhPlZR/+/qUX08kLn6aISgNIulQAUu9LQINl9vIQ4E3X/WCMhBGPDdew4WqwmiVHapBKA0jn+/vrFxhw8e4HSylqTjP9v5O/hj7+HdbXF2HAdReXVHM5XcwUrV6aGgSqds3q11hF0yq5L4/91Jtvn6lldb4vGSqAxCTFd3hfoEwDA/44UEhGkbalrxXqpMwDFriQfK6afrgdB3m5ah3KZrpaDbi7Ex53+uh7qOoDSJpUAlM554QX9YkPq6hv48XgJo008+ueFpBd4Icn62mLcQB27j5dQU6dmCVNapxKA0jnJyfrFhhzIvUBZdZ1J+/8BknOSSc6xvrYYG66joqae1NPntQ5FsVIqASh2o7H+z+j+fhpHYhmj+/vjIFB3BStXpBKAYjeSjxUTHuhJoJd19f+bi7e7M8NCfdR1AOWKtJwUPkwIsUUIcUgIcUAI8ZRWsSjdX219A3tOlHTbu3+vZFy4jtTT5ymrUuWhlZa0PAOoA56RUg4BRgOPCyGGahiP0hGhofrFRuzPKaWipt7k/f8AoT1DCe1pnW0xNlxHfYNkd3aJ1qEoVkjLSeHPAmcNj8uEEIeAEOCgVjEpHfDvf2sdQYfsutT/b/oE8O87u94WpioH3dzVV/ng5uzAjqNFTBzayyzHUGyX0QlACOEI9Gr6HinlKVMEIYToC4wAdptif4rSXPKxYgYHeeHXw0XrUCzK1cmRa/v5W9WF4NraWnJycqiqqtI6lG7Hzc2N0NBQnJ2djdreqAQghPgN8AcgH2gcVCyBYZ0Jstm+PYHVwDwp5YVW1s8B5gD06dOnq4dTTGXePP3fxYu1jMIo1XX1pJws4e6R5vn8zNs4D4DFkxabZf9ddV24jje/OUReaZVV3ACXk5ODl5cXffv27TbzMVsDKSXFxcXk5OTQr18/o95j7BnAU0CElLK409G1QgjhjP7Lf7mU8svWtpFSLgOWAcTGxqrKVtYiNVXrCIyWdrqUqtoGs/T/A6TmpZplv6bSWBbih6NF3HWN9tcqqqqq1Je/GQgh8Pf3p7Cw0Oj3GHsR+DRQ2qmorkDo/+9/CBySUv7FlPtWlKaSjxUjBIzuZ18jgBoNDvLCv4d1lYdWX/7m0dF2NfYMIBvYKoT4L1Dd+GIXv7jHAvcD6UKIVMNr86WU33Rhn4rSQnJ2EUODe+LtYVy/aHfj4CCIC9ex42gRUkr15atcYmwCOGVYXAxLl0kpdwDqk6iYVVVtPT+dOs8vR1+ldSiaGhfuz7q0XI4WlDOwl6oOqugZlQCklK8BCCG89E+lKjJu7wYN0joCo/x06hw1debr/wcY5N/1tjB1OejmGq8D7DhapBKAjaivr8fR0fGKz03BqGsAQogoIcQ+IAM4IITYK4SINGkkim1Ztky/WLldx4pxEDCyn/nq/yybuoxlU7vWFqYuB91cqK8Hff09rOo6gJZOnDjB4MGDefjhh4mKiuLee+8lKSmJsWPHMnDgQH788Ud+/PFH4uLiGDFiBHFxcWRlZQHwl7/8hYceegiA9PR0oqKiqKioaPU45eXlPPjgg0RHRzNs2DBWG+bR+Pzzz4mOjiYqKornnnvu0vaenp688sorjBo1iuTk5BbPTc3YLqBlwG+llFsAhBDxwP8BcSaPSFFMKDm7mOgQb3q62Wf/f1Njw3V8nZpLbX0Dzo7WUwYsPjG+xWuzImcxd+RcKmormLx8cov1CTEJJMQkUFRRxIyVMy5btzVhq1HHPXr0KP/5z39YtmwZI0eO5LPPPmPHjh2sXbuWhQsX8q9//Yvt27fj5OREUlIS8+fPZ/Xq1cybN4/4+HjWrFnDm2++ydKlS/HwaH1+6QULFuDt7U16ejoA586dIzc3l+eee469e/fi6+vLzTffzFdffcUdd9zBxYsXiYqK4vXXXwdo8dzUjP0U9Gj88geQUm4FepglIsU2zJmjX6xYpaEU8mgzz/87Z90c5qyz7rYAuG6gjvLqOtJUeWgA+vXrR3R0NA4ODkRGRnLjjTcihCA6OpoTJ05QWlrKzJkziYqK4umnn+bAgQMAODg4kJiYyP3338/48eMZO3bsFY+RlJTE448/fum5r68ve/bsIT4+noCAAJycnLj33nvZvn07AI6Ojtx1112Xtm/+3NSMHgUkhHgZ+NTw/D7guHlCUmzC4cNaR9CulJMl1NZLsxeAO1xs/W0B+mkwhaE8dGxf6ymJ3dYvdg9njzbX6zx0Rv/ib87V1fXSYwcHh0vPHRwcqKur4+WXX2bChAmsWbOGEydOEB8ff2n7I0eO4OnpSW5ubpvHaG3UlZRXvp3Jzc3tsn7+5s9NzdgzgIeAAOBLYI3h8YPmCkpRTCH5WDFODoKRVvRlpyVvD2eGhXir6wBGKi0tJSQkBIDExMTLXn/qqafYvn07xcXFrFq16or7uPnmm3nvvfcuPT937hyjRo1i27ZtFBUVUV9fz+eff8748ePN9t/RFqMSgJTynJTySSnl1VLKEVLKp6SU58wdnKJ0RXJ2McNCvenhqlnNQ6szNlzHvlPnKa+u0zoUq/f73/+eF154gbFjx1JfX3/p9aeffpq5c+cyaNAgPvzwQ55//nkKCgpa3cdLL73EuXPniIqKYvjw4WzZsoXg4GDeeustJkyYwPDhw7n66qu5/fbbLfWfdZk2/2UIIRZLKecJIdahr/1zGSnlNLNFpihdUF5dx/6cUh4b31/rUKzKuHAdH2w9xo/Hi7lhsP1WB+3bty8ZGRmXnjf9hd903eEmXZ0LFiwA4KOPPrr0WlhYGEePHr3icTw9Pfnkk09avH7PPfdwzz33tHi9vLy8zeem1t5Po8Y+/0VmjUKxPTExWkfQpj0nSqhvkIzprzP7sWKCYrq8D3OVg27u6qt8cXVyYMcR+04Ail6bCUBKudfwMEZK+W7TdYYZvLaZKzDFyll5FdBdx4pxdhRcc5Wv2Y9lrVVAW+Pm7Mi1/fzUdQAT+/jjj3n33cu+Ihk7dizvv/++RhEZx9jO0QeAd5u9ltDKa4piFZKzixkR5ou7i/lGUNiqseE6/rghk4ILVQT21L48dHfw4IMP8uCDtjcups2LwEKIXxj6//sJIdY2WbYAJi0NrdiY++7TL1boQlUtGWdKzT7+v9F9X97HfV9aZ1u0Zlxjeehj6izA3rV3BrAT/bSNOuCdJq+XAfvNFZRiA3JytI7gin7MLqFBYrEJ4HMuWG9btGZocE98PZzZcaSY6SO0nx9A0U571wBOAieBMZYJR1G6Ljm7GBcnB0b08dE6FKvUWB76B1Ue2u4ZWwxutBBijxCiXAhRI4SoF0K0mL5RUaxB8rFirunji5uz6v+/knHhOvIuVHGs8KLWoSgaMvZO4PeAXwBHAHfgYeDv5gpKUTrrfEUNh/IumLX8szmYuxpoc+OaTBOp2C+jSwJKKY8CjlLKeinlx8AE84WlWL0xY/SLldmVXYKUWDQBjAkdw5jQrrWFpRNAmJ8Hffw82KESwBWtXbuWP/7xj1qHYVbGDgOtEEK4AKlCiLfRXxhW1UDt2VtvaR1Bq5KPFeHu7MjwUB+LHfOtidbZFu0ZG65jfVoudfUNOFlReWhrMW3aNKZNs2yxg7q6OpycnK743NSM/b9+P+AIPAFcBMKALtcoFUJ8JIQoEEJktL+1orQvObuY2L6+uDipL7T2jAvXUVZdR1pOqbaBxMe3XD74QL+uoqL19Y2lG4qKWq4zgjETwiQmJvLEE08AkJCQwJNPPklcXBz9+/dvswAcwNtvv010dDTDhw/n+eefByA1NZXRo0czbNgwpk+fzrlz5wz/+fHMnz+f8ePH8+6777Z4bk7GFoM7KaWslFJekFK+JqX8raFLqKsSgUkm2I9iaXfdpV+sSFF5NYfzyxltoeGfje5aeRd3rbSutjBG3AB/hLDf6wBHjx7lqaeeYv/+/WRmZl6aEGbRokUsXLiwxfZnz55lx44drF+//tKXems2bNjAV199xe7du0lLS+P3v/89AL/85S/505/+xP79+4mOjua111679J7z58+zbds2nnnmmVafm0t7xeDSaaUIXCMp5bCuHFxKuV0I0bcr+1A0Umx99wHuytbHZOkLwMUV1tcWxvDt4UJUb292HC3iyRsHahfI1q1XXufh0fZ6na7t9W1onBAGaHVCmObuuOMOHBwcGDp0KPn5+Vfcb1JSEg8++OClWcL8/PwoLS3l/Pnzl8o+P/DAA8ycOfPSe2bPnn3ZPpo/N5f2OpemWCQKRTGB5GPF9HBxJDrEW+tQbMbYcB0f7sjmYnWd3ZXNbm9CmLa2b2tSl87cW9GjR482n5tLm11Ahq6fk4YbwgAGGh4XACVmjw4QQswRQqQIIVIKCwstcUjFRiVnFzOyn59VzXdr7caF66itl/x4wiL/nO3CzTffzEcffXRpoviSkhK8vb3x9fXlf//7HwCffvqpZpPANGVUyhdCPALMAfyAAUAo8A/gRvOFpielXIZ+UnpiY2OvnHYVu5Z/oYrswovMjg3TOpROsVQ56OYaL5j/cKSICRGBmsTQ3UyaNInU1FRiY2NxcXFh8uTJLFy4kE8++YTHHnuMiooK+vfvz8cff6x1qIi2TmUubSREKnAtsFtKOcLwWrqUMrrLAeivAayXUka1t21sbKxMSUnp6iEVUzBMjsHLL2sbh8HXqWd46otU1j4xlmEWHAIKsGCbvi1eHm8dbdFR9/5zF8XlNWycd71Fjnfo0CGGDBlikWPZo9baVwixV0oZ23xbYzv9qqWUNY39WkIIJ9q4OGwsIcTnQDygE0LkAH+QUn7Y1f0qFmAlX/yNko8V4+XmRGRvy/f/2+oXf6Ox4Tre3phFYVk1AV6u7b9B6TaMTQDbhBDzAXchxE3AXGBdVw8upfxFV/ehKKDv/x/Vzw9HB1XYrKOuCw/gbbLYeayI22NCtA7HZqSnp3P//fdf9pqrqyu7d+/WKKKOMzYBPIe+/k868CjwDfBPcwWl2IBbb9X/3bBB2ziA3POVnCyu4P7RV2ly/FuX69tiw73at0VnDO3dEx8PZ3YcUQmgI6Kjo0lNTdU6jC5pNwEIIRyA/YY++v8zf0iKTais1DqCS5KPaTP+v1FlrfW0RWc4OgjiBvir8tB2qN3xclLKBiBNCNHHAvEoSoclZxfj4+HMkKCeWodis8aG68gtreJ4kSoPbU+M7QIKBg4IIX5EXwsIACmlZSslKUorko8VM7qfPw423P/fWAk0JiFGk+M3LQ/dP8BTkxgUyzP2jpnX0N8V/Dr6qSEbF0XR1OmSCs6cr7S5+v/NWbocdHN9/DwI9XVX5aFb0bQoXEclJCS0WzguMTGR3NzcDu33xIkTREW1O3K+XcZeA3jfmHH6ih2ZYh1VQrTu/weYMsg62qIrhBCMC9fx3/Sz1DdINZrKghITE4mKiqJ3794WP3a7CUBK2SCESBNC9JFSnrJEUIoNePZZrSMA9P3/Ok8XBgZq123xbJx1tEVXjQ3X8cWe06SfKSUmzMdix02MT7ziupiEmDa7xVp7r7F3Vf/rX/9i0aJFCCEYNmwYs2bN4o033qCmpgZ/f3+WL19Or169LntPfn4+jz32GNnZ2QAsWbKE3r17M2XKFDIy9FXtFy1aRHl5Oa+++upl73399ddZt24dlZWVxMXFsXTpUlavXk1KSgr33nsv7u7uJCcnc/DgQX77299SXl6OTqcjMTGR4OBg9u7dy0MPPYSHhwfjxo0z6r+xPcZ2ATVeA9gshFjbuJgkAkXpJCklyceKGdXfX41cMYGxdjRN5IEDB3jzzTf5/vvvSUtL491332XcuHHs2rWLffv2cffdd/P222+3eN+TTz7J+PHjSUtL46effiIyMtLoYz7xxBPs2bOHjIwMKisrWb9+PTNmzCA2Npbly5eTmpqKk5MTv/nNb1i1atWlL/wXX3wRgAcffJC//e1vJCcnm6wdjL0I/Fr7myh2pXHijU6W4jWFE8UV5F2oYoyF6/83F58YD8DWhK2axtFVfj1ciOzdk+2HC3l8QrjFjtuVOkidfe/333/PjBkz0On0Sc/Pz4/09HRmz57N2bNnqampoV+/fq2+71//+hcAjo6OeHt7X5rYpT1btmzh7bffpqKigpKSEiIjI5k6depl22RlZZGRkcFNN90EQH19PcHBwS3KSd9///1sMME9OMZOCLMNyAS8DMshw2uKopnvMwsAuG6gTuNIuo/4iABSTp7j3MUarUMxq9bud/jNb37DE088QXp6OkuXLqWqqsqofTk5OdHQ0HDpeWvvq6qqYu7cuaxatYr09HQeeeSRVreTUhIZGUlqaiqpqamkp6fz3Xffme3+DKMSgBBiFvAjMBOYBewWQswweTSK0gFJB/MZ1MuTq/zV9NSmcmtUMPUNkk2HrjzhSXdw4403snLlSooNExuVlJRQWlpKSIj+TuhPPvnkiu9bsmQJoP91fuHCBXr16kVBQQHFxcVUV1ezfv36Fu9r/LLX6XSUl5dfNjLIy8uLsrIyACIiIigsLLzUzVNbW8uBAwfw8fHB29ubHTt2ALB8+XJTNIPRXUAvAiOllAUAQogAIAloe3yTophJaUUtP54o4dHr+2sdikloVQ66ucjePQnxcefbjDxm2WhpbWNERkby4osvMn78eBwdHRkxYgSvvvoqM2fOJCQkhNGjR3P8+PEW73v33XeZM2cOH374IY6OjixZsoQxY8bwyiuvMGrUKPr168fgwYNbvM/Hx4dHHnmE6Oho+vbty8iRIy+tS0hI4LHHHrt0EXjVqlU8+eSTlJaWUldXx7x584iMjOTjjz++dBH4lltuMUk7GFsO+rLSz4ahoWmmKAfdEaoctBXR+BpAY/nnL+fGcXUfX01iaNRdrgE0WrD+IJ8mn+SnV27C0wyzhKly0OZljnLQG4UQ3wKfG57PRl8QTrFXs2ZpevikQwXoPF2IsXDt/9bMitS2LUzt1qggPtxxnO8zC5g23PJj0xXLaW9S+HCgl5Tyd0KIO4FxgACSAdN0Qim2ae5czQ5dU9fA1qwCbo0KsoryD3NHatcW5nB1H18CvFz5NiNPJYBurr2LwIuBMgAp5ZdSyt9KKZ9G/+t/sXlDU6xaRYV+0cCeEyWUVdUxcUiv9je2gIraCipqtWkLc3BwENw8tBdbsgqoqq03yzGM6XpWOq6j7dpeF1BfKeX+Vg6SYpjK0Sb8cUMm69JyCfZ2I8jbjd4+7gT1dCPY241gH3eCvd3Qebqq2987YvJk/V8NrgFsOpiPq5MD46xk+Ofk5fq26C7XAEA/Gmj57lNsP1zIzZFBJt23m5sbxcXF+PurG/hMSUpJcXExbm5uRr+nvQTQ1p7cjT6KxoYEe5F/wY+zpZVknCnlu4P51NQ1XLaNk4OgV099ggg2LEHe7vRukjQCPF2tosvBnkkpSTqUz7hwHR4upr9AqeiN6u+Ht7szGw/kmTwBhIaGkpOTQ2FhoUn3q+iTa2hoqNHbt/cvaI8Q4hEp5WUTwQghfgXs7UR8lxFCTALeBRyBf0op/9jVfbbm9piQy2Y6klJyrqKW3POV5JVWcfZCFWcbH5dWkXGmlE0H86luliQGB3mR+OC1BHkbn2EV08rKLyPnXKVF71S1BK3LQTfn7OjAxCG92HQwj5q6BlycjK0aY8S+nZ1bvctWsbz2EsA8YI0Q4l5+/sKPBVyA6V05sBDCEXgfuAnIQZ9s1kopD3Zlv0YeG78eLvj1cCEqpPVJxBuTxNlSfWI4UVzBXzcdZsY/drL84VHq5iONJB3U36B04+BAjSMxLWtLAKAfDbT6pxx2ZRdz/aAArcNRzKDNtC6lzJdSxqGvBXTCsLwmpRwjpczr4rGvBY5KKbOllDXAF8DtXdynyTQmicje3tw4pBe/GtePzx4ZRXl1HTP/kUxWXpnWIdqlpEMFDA/zIbCnOgszt3EDdXi4OLLxQFf/qStdVVBmXFmKjjK2FtAWKeXfDcv3Jjp2CHC6yfMcw2tWa1ioDysfHQPA7GXJpJ4+r21AWkpI0C8WVFBWRerp89w0xLp+/SfEJJAQk6B1GCbn5uzIhMGBfHcgj/oGNWpHC1JKPt11kuv+tIVth01/zcR0HXsd19rV1BafMiHEHCFEihAixRouGg3q5cWqx+LwcnPi3v/bdWlCErujQQL4/pC++NuNVjL8s1F3TQCg7wYqKq9h70njKl4qpnOxuo55K1J5+asMxgzwZ9gVuqu7QssEkAM0LTYSCrSYF01KuUxKGSuljA0IsI5+yD7+Hqx6LI7ePu488PGPbO7mhbNaVVSkXywo6VA+IT7uDA7ysuhx21NUUURRRfesoR8fEYiLkwMbM1Q3kCVl5ZUx7b0drEvL5Xe3RPDRAyPx7eFi8uNomQD2AAOFEP2EEC7A3YDNTDLTq6cbKx4dw+AgLx79dC9fp57ROiTLmjFDv1hIZU09/ztSxE1De1nd2PEZK2cwY2X3LI7r6erE9QN1fHsgT928ZSGr9+Zw+/s7KK2sY/nDo3l8QrjZhp9rlgCklHXAE8C3wCFgpZTygFbxdIZfDxeWPzyKa67yZd6KVJbvPql1SN3WjqNFVNc1WM3dv/ZkUlQwZ85Xkn6mVOtQurWq2nqeW7WfZ/6TRkyYD988Nc7sc11reieNlPIbbLyonJebM588dC1zl//Ei2syKKuq47HxA7QOq9tJOpiPl6sT1/bz0zoUs7CWctCtmTgkEEcHwcaMPIZZQfG97uh40UV+/e+9ZOaV8cSEcOZNHIiTo/l/n2vZBdRtuDk78o/7rmHKsGD+uCGTP3+bqU6XTaihQbI5s4DxEQEmvSFJMY6Phwtj+vuzMUN1A5nDf/efZerfd5B/oYqPHxzJs7dEWOTLH1QCMBkXJwfevXsEv7g2jPe3HOMPaw/QoIbOmURaznmKyqu5aajq/tHKpKggsosucqSgXOtQuo2augZeXXuAxz/7iYG9PPnvk9cxIcKyQ5xVMRUTcnQQLJwejZebM8u2Z1NWVcefZwyzWDa3qF//2mKHSjqUj6ODIH6QdY3/b/TrWMu1hVZuHtqLl7/OYEN6HoN6WdcoLFuUc66Cxz/bR9rp8zw0th/P3zpYk7NblQBMTAjBC7cOpqebE4u+O0x5dR1//8UI3JwdtQ7NtGbPttihkg4WMLKvL94ezhY7ZkfMjrJcW2glsKcb1/TxZeOBPJ6aOFDrcGza95n5PL0ijYYGyT/uu5pJUcGaxdINf5pqTwjBEzcM5LVpkWw6mM+vPtnDxeo6rcMyrdOn9YuZnSquICu/zKpH/5wuPc3pUvO3hdYmRQVx6OwFThZf1DoUm1RX38CfNmbyUGIKob7urH9ynKZf/qASgFk9ENeXd2YOJ/lYMfd9uJvSilqtQzKd++/XL2aWZLjJzpr7/+9fcz/3rzF/W2jtFkNZaHVTWMflX6jinn/uZsnWY/zi2j6s/nWcVRSUVAnAzO66JpQP7r2GA2cuMHtZstmKOnVXSYfyGRjoaRX/WMwpNTH1UkVQaxXm50F0iLcqDtdBPxwt4ra//Y/0nFL+Ons4b90ZbTVdwioBWMCkqCA+TIjlZHEFD3+SQl19Q/tvUiitrOXH4yVMtOJf/6ZiCwkA9J/lfafOk1eqfsgY48fjJfzyox/x8XBh7RNjmT7C+MlaLEElAAu5bmAAb88Yxv6cUj5JVncMG2Pb4ULqGqRV9//bm8ZuoG/VWUC7SitreXpFKqG+7qyZG8dAKxw9pRKABU0ZFkx8RADvfJfFmfOVWodj9ZIO5uPfw4WYMB+tQ1EMwgM9GRjoqa4DGOGVrzPIu1DF4tkxeLlZ5wg2lQAsSAjBgtujkBL+8HWGbd9V+cwz+sVMausb2JJVwA2D9WUIrNkzY57hmTHmawtrMykqiN3Hiym5WKN1KFbrq31n+Do1l6duHMiIPr5ah3NFKgFYWJifB0/fNJCkQwW2/Stq6lT9YiZ7jpdQVlVnE/3/UyOmMjXCfG1hbW6JDKJBwqaDNvz5NaPTJRW8/FUGsVf5MjfeuuuCqQSggYfG9mNocE/+sPYAF6psdGhoVpZ+MZNNh/JxcXLguoE6sx3DVLKKssgqMl9bWJvI3j0J83O37R8wZlJX38DTK1IB+OvsGKuvAmDd0XVTTo4OvHVnNEXl1by9MVPrcDrn0Uf1ixlIKUk6lM+4cB0eLtZ/s/qj6x/l0fXmaQtrJIRgUmQQO44W2e4PGDNZsvUYKSfPseCOKML8PLQOp10qAWhkeJgPD8T1ZfnuU2q6vWYO55dzuqTSrkb/JGxNsOqS0M1Nigqitl6yJbNA61CsRurp8yzefIRpw3tzxwirnt78EpUANPTMzREE9XRj/pfp1Kp7Ay5pvPv3Riub/F352YgwXwK9XFU3kMHF6jqe+mIfQT3dWHBHlNbhGE0lAA15ujrx+u1RZOWXsWx7ttbhWI2kQ/kMD/WmV083rUNRrsDBQXBLZBBbswqprKnXOhzNvbbuAKdLKvjr7Bi83a1zyGdrVALQ2E1DezEpMoi/bT7CiSJVZKugrIrU0+e50Y66f2zVpKggKmvr2Xa4UOtQNLUh/SwrU3L4dfwAm5uxTpMEIISYKYQ4IIRoEELEahGDNXl1WiTOjg68+FW67dwb8NJL+sXEtmQWICU21f//0vUv8dL1pm8Lazeqnx8+Hs52fVfw2dJKnv8ynWGh3sybOEjrcDpMqzOADOBOYLtGx7cqQd5uPDcpgh+OFrNm3xmtwzHOxIn6xcQ2HSwgxMedIcHWd9v8lUzsP5GJ/U3fFtbOydGBm4b0IulQPjV19ncNq6FB8szKNGrqGnj37hE4W/mQz9ZoErGU8pCU0n4GThvh3lFXMaKPD2/895Bt3GGZmqpfTKiypp4dRwuZOCQQIaz77t+mUvNSSc1L1ToMTUyKCqKsqo6dx4q0DsXi/rkjm53HivnD1KH009lmtVrbS1ndlIOD4K07o7lQWcvCbw5pHU775s3TLyb0w9EiqmobbOLu36bmbZzHvI3zurQPW6kG2tzYcB2erk52NxroQG4pf/42i1siezF7ZJjW4XSa2RKAECJJCJHRynJ7B/czRwiRIoRIKSzs3hebBgf15JHr+7Nqb45d/qJKOpSPp6sTo/r5ax2KxdlqAnBzdmTC4EC+O5hPfYONXL/qosqaep76IhW/Hi788c5hNnW22pzZEoCUcqKUMqqV5esO7meZlDJWShkbEBBgrnCtxlM3DqSPnwcvrsmgqtZ+htc1NEg2ZxYwPiJAk8mxlc6bFBlEycUa9pwo0ToUi1j4zSGOFpSzaOZwfHu4aB1Ol6h/aVbGzdmRN6dHcbzoIu9vOap1OBaz/0wphWXVTFQ3f9mc+IgAXJ0c7KIbaPOhfD7ddZKHx/XjuoG2/4NUq2Gg04UQOcAY4L9CiG+1iMNaXTcwgOkjQvjHtmMcyS/TOhyLSDqYj6ODYEKESgC2poerE9cPCmBjRh4N3bgbqLCsmt+v2s/gIC9+NylC63BMQpNKW1LKNcAaLY5tK166bQhbsgp44ct0Vj46Bgdrq4m/cKFJd5d0KJ/Yq3zx8bC9U+qFN5q2LWzRpMggNh3MZ/+Z0m45gY+Ukt+tSqO8uo7P54zG1ck65vTtKtUFZKX8PV2ZP3kIKSfP8cWe01qH01JcnH4xgdMlFWTmlXGTjY3+aRQXFkdcmGnawlZNHNILJwfBhoyzWodiFv9KPsnWrELmTx7CICuc2rGzVAKwYjOvCWV0fz/e2nCIgjIrm4R75079YgI/F3+zzQSw8/ROdp42TVvYKm8PZ8YM8OfbjDzbuZvdSIfzy1j4zSEmRATwyzFXaR2OSakEYMWEELw5PZrq2gZeX3dQ63AuN3++fjGBpEP5hAd62uzNNPM3z2f+5q61ha2Vg27NrVHBnCiuIKsbXbeqrqvnyc/34enqxNszhtv0kM/WqARg5QYEePL4hHDW7z/bLWuvX6iqZXd2iU3V/lFad9PQXggBG9K7z2igP2/MIjOvjLdnDCPAy1XrcExOJQAb8Fh8f8IDPXnpqwwqauq0DsektmUVUtcguWmoGv1j6wK8XBl5lV+3KQ63JbOAf+44zn2j+9hs92R7VAKwAa5Ojrx1ZzRnzlfy102HtQ7HpJIO5ePXw4WYMF+tQ1FMYFJUEJl5ZRy38dLmZ0sr+e3KVAYHefHSbUO1DsdsVAKwESP7+vGLa8P46IcTZJwp1Tock6itb2BLZgE3DA7E0dqGuSqdcktUEIBN3xRWV9/AU5+nUl3XwPv3Xo2bc/cY8tkalQBsyPOThuDr4cLvVu3Xvvzu4sX6pQv2nCjhQlWdzff/L560mMWTFmsdhlUI8XFneJgP/951kjIbnTB+cdIRfjxRwpvToxgQ4Kl1OGalEoAN8fZwZuH0KA6dvcC7mzXuCoqJ0S9dkHSwABcnB64bqDNJSFqJCYohJihG6zCsxsu3DeFsaSVvrLeBqrbN/O9IIe9vPcqs2FCmjwjVOhyzUwnAxtwcGcTMa0JZsvUYe09qWHwrKUm/dJKUkk2H8hg7wJ8erprckG4ySdlJJGV3vi3AdquBtia2rx+Pjh/AipTTbDqYr3U4Riu4UMW8L1IZGOjJa9NsZ2L3rlAJwAa9MnUowd7u/HZlGherNRoV9MYb+qWTth0u5HRJJTcNDTJhUNp4Y/sbvLG9820B3SsBADw9cRBDgnvywpf7KS6v1jqcdtU3SJ76IpWKmnrev+dq3F26b79/UyoB2CAvN2femTWcUyUVvLXB9k6zL1bX8eKaDPoH9ODOq0O0DkcxAxcnB/4yazgXKuuYv8b657r+2+YjJGcX8/rtkQzsRqUe2qMSgI0a3d+fX43tx793nWJrlm3dIPaXTYc5c76SP945rFuPsLB3Q4J78tubB/HtgXy+/Ml657reebSIv31/hDuvDmFmrO3O7tUZKgHYsGdviWBgoCe/X7Wf8xU2MI8wkHb6PB//cJx7R/Xh2n5+WoejmNkj1/VnZF9fXl17gDPnK7UOp4XCsmqeWpFKf10PFtxuH/3+TakEYMPcnB356+wYSi7W8PLXB7QOp1219Q08t3o/AV6uPHfrYK3DUSzA0UHwzswY6qXk2ZVpVjVfQH2D5OkVqVyorOX9e6+2+cEInaESgI2LCvFm3sSBrEvLZW1aruUOvHSpfumAZduzycwrY8HtUfR0czZTYJa3dMpSlk7pWFvYkz7+Hrw8ZSjJ2cUk7jyhdTiXfLDlKDuOFvHatEgGB/XUOhxNqATQDTw2fgAxYT68/FUGeaUWKhsdEaFfjJRdWM67m48wOTqImyNtf+RPUxG6CCJ03WOGKHO5e2QYNwwO5E8bMzlaoH210F3Zxfw16TC3x/Rm9kj76vdvSiWAbsDJUT/iorqunt+v3m+ZERfr1ukXIzQ0SF74Mh03JwdenRZp5sAsb13WOtZlGdcWV9IdykG3RQjBH++KxsPFkadXpFFbr92d7MXl1Tz1xT6u8u/Bm9Oju12J547Qak7gPwshMoUQ+4UQa4QQPlrE0Z30D/DkxclD2H64kOW7T5n/gO+8o1+MsCLlNLuPlzB/8hACvdzMHJjlvZP8Du8kG9cW9izQy42F06NJP1PKe98f1SSGhgbJb1emca6ilvfuGYGnHfb7N6XVGcAmIEpKOQw4DLygURzdyn2jr+K6gTre/O8hq6nGmH+hioXfHGJ0fz+7PtVW9G6NDmb6iBDe23KUtNPnLX78pduz2Xa4kJenDCWyt7fFj29tNEkAUsrvpJSNt7DuArp/0Q0LEELw5xnDcXYUPLMylToNT7Mb/eHrA1TXNfDWncPs+lRb+dmr0yIJ9HLl6ZWpVNbUW+y4KSdKWPRdFrdFB3PfqD4WO641s4ZrAA8BG7QOorsI8nZjwR1R/HTqPEu3Z2say8aMPDYeyGPexIE2O92jYnre7s78ecZwsgsv8qeNmRY55rmLNfzm832E+Ljz1l323e/flNkSgBAiSQiR0cpye5NtXgTqgOVt7GeOECJFCJFSWFhornC7lWnDe3PbsGAWJx3mQK42cweUVtbyytcZDAnuySPX9dckBsV6jRuoIyGuL4k7T7DjSJFZjyWl5Nn/pFFcXsP791zdrYYgd5XQqkaHEOIB4DHgRillhTHviY2NlSkpKeYNrJs4d7GGWxZvx9fDha+fGGv6kgunT+v/hrXerz9/TTpf/HiKrx4fy7BQH9Me28qcLtW3RZi3usbREZU19dz29/9RWVPPxnnX4+1uni/m/9uezZvfHOLVqUNJGNvPLMewdkKIvVLK2OavazUKaBLwHDDN2C9/pWN8e7jwpxnDyMovM880kmFhV/zy351dzGe7T/Grcf26/Zc/6L/4u/rl392qgRrD3cWRv8yKoaCsmtfWmudO9p9OneNPGzOZFBnEA3F9zXIMW6bVNYD3AC9gkxAiVQjxD43i6NYmRARyz6g+LPtfNruzi0278xUr9EszVbX1vPBlOqG+7jx90yDTHtNKrchYwYqMlm3REfaYAABiwnx4fEI4X+47w4b0sybdd2lFLb/5bB9B3m78aYYahNAarUYBhUspw6SUMYblMS3isAcvTh5CHz8PnvlPGuWmnDtgyRL90sx73x8lu+giC6dH4+FiH2Osl6QsYUlKy7ZQjPObG8KJDvFm/pp0Csq6fif7xeo6Nh3MZ86nKRSUVfHePVebrXvJ1lnDKCDFjHq4OvHOzOHknq/kjfUHzXqszLwL/GPbMe68OoTrBwWY9VhK9+Hs6MBfZw/nYk09L6zu3NwBx4su8tGO49z/4W5GvL6JR/6VwoHcCyy4PYqYMB/TB91N2MdPNDvXOEXfkq3HmDikFxOHmn4S9voGyXOr0/F2d+bl24aafP9K9xYe6MVzkwazYP1BVuw5zd3Xtj1Ov6q2nh+Pl7Alq4AtmQWcKK4w7MeTB+KuYkJEILF9/XBxUr9x26ISgJ14euIgtmYV8vyX+/m2z/X4e7qadP+f7DxB2unzvHt3DL49XEy6b8U+PBjXl6SD+SxYf5C4ATr6+Htctv7M+Uq2ZhWwJbOQH44WUVlbj6uTA3ED/HloXD/iBwW2eI/SNpUA7ISLk/40e9rff+DFNRksue9qk10UO11SwaLvspgQEcC04b1Nsk/F/jg4CBbNGs6kv27n2f+k8enD15J66jxbsgrZkllAVr6+imiorzszY0OZEBHI6P7+djN/rzmoBGBHBgf15JmbB/HWhkyeW72fydHBjO7v37l7BFatAvQ32bz4VQYAb9hpZcVVs1ZpHUK3EeLjzqvTInnmP2kMe/U7qusacHIQXNvPjxevGcKEwQEMCPC0y8+ZOagEYGcevq4/xwrLWZuWy8qUHNycHRg7QEf84EAmRAQQ6mvkKbROB8DX+86w/XAhr04dSoiPuxkjt146D12X99GdS0F31J1Xh3CkoJxzF2uYMDiAseE6vNTdu2ah2Z3AnaHuBDadqtp6dh8vYUtmAd9nFnCqRH8RbWCgJzcMDiQ+IpDYvr44O17hIlpiIuXVdVxfcBV9/DxY/es4HB3s81dZYmoiAAkxCZrGoShXcqU7gVUCUJBSkl10kS2ZBWzNKmT38WJq6yVerk6MG6hjwuBA4gcFENizSS3/+HiOFpQz6Y7X+e+T1xER5KXdf4DG4hPjAdiasFXTOBTlSq6UAFQXkIIQggEBngwI8OTh6/pTXl3HD0eLLo242JCRB0BUSE9uiAgkfnAg/SpqKSqvZm78ALv+8lcUW2ZbCSArC+LjL39t1iyYOxcqKmDy5JbvSUjQL0VFMGNGy/W//jXMnq0vbnb//S3XP/MMTJ2qP/ajj7Zc/9JLMHEipKbCvHkt1y9cCHFxsHMnzJ/fcv3ixRATA0lJ8MYbLdcvXaqfe3fdutZn4Pr0U31NnhUrWr0zl1Wr9P31iYn6pblvvgEPD/jgA1i5EgBP4BbDIrds4dDZMs69vhCfL7+jrKqWGsCxIJuhDo6MuCFcv58FC2Dz5sv37e8Pq1frH7/wAiQnX74+NBT+/W/943nz9G3Y1KBBsGyZ/vGcOXC4WU2jmBh9+wHcdx/k5Fy+fswYeOst/eO77oLiZuUwbrwRXn5Z//jWW6Gy8vL1U6bAs8/qHzf/3MGlz55rdT1/+ms6JDbbRn32OvzZu8zWrfq/ixbB+vWXr3N3hw2GKvJ2/Nnr9PeegW0lAMXihBAM7d0TwnWQ2ZO6Bsn5ilqcix1wdnbEyUkNwVMUW6WuASid0/irpPFXmh1T1wAUa6euASim9c03WkdgNb65t+tt0VgJNCYhpsv7UhRjqQSgdI6HuuW+kYdz19tCJQBFC6pSktI5H3ygXxQ+2PMBH+xRbaHYHpUAlM5ZubL1kRt2aOWBlaw8oNpCsT0qASiKotgplQAURVHslEoAiqIodkolAEVRFDtlUzeCCSEKgZNax9EOHVCkdRBGUHGalq3ECbYTq4rTdK6SUraYqNumEoAtEEKktHbHnbVRcZqWrcQJthOritP8VBeQoiiKnVIJQFEUxU6pBGB6y7QOwEgqTtOylTjBdmJVcZqZugagKIpip9QZgKIoip1SCcAEhBAzhRAHhBANQojYJq/3FUJUCiFSDcs/rDFOw7oXhBBHhRBZQohbtIqxNUKIV4UQZ5q0YytTIGlHCDHJ0G5HhRDPax3PlQghTggh0g1taFUTawghPhJCFAghMpq85ieE2CSEOGL466tljIaYWovTqj+fbVEJwDQygDuB7a2sOyaljDEsj1k4ruZajVMIMRS4G4gEJgEfCCGsbaqvvzZpR6uZjMDQTu8DtwJDgV8Y2tNaTTC0obUNW0xE/9lr6nlgs5RyILDZ8FxribSME6z089kelQBMQEp5SEqZpXUc7WkjztuBL6SU1VLK48BR4FrLRmezrgWOSimzpZQ1wBfo21PpACnldqCk2cu3A58YHn8C3GHJmFpzhThtlkoA5tdPCLFPCLFNCHGd1sFcQQhwusnzHMNr1uQJIcR+wym45l0BTdhC2zWSwHdCiL1CiDlaB2OEXlLKswCGv4Eax9MWa/18tkklACMJIZKEEBmtLG392jsL9JFSjgB+C3wmhOhphXGKVl6z6PCwduJeAgwAYtC36TuWjK0dmrddB4yVUl6NvrvqcSHE9VoH1E1Y8+ezTWpKSCNJKSd24j3VQLXh8V4hxDFgEGC2C3CdiRP9r9awJs9DgVzTRGQcY+MWQvwfsN7M4XSE5m1nLCllruFvgRBiDfruq9auW1mLfCFEsJTyrBAiGCjQOqDWSCnzGx9b4eezTeoMwIyEEAGNF1OFEP2BgUC2tlG1ai1wtxDCVQjRD32cP2oc0yWGf/yNpqO/mG0t9gADhRD9hBAu6C+mr9U4phaEED2EEF6Nj4Gbsa52bM1a4AHD4weArzWM5Yqs/PPZJnUGYAJCiOnA34EA4L9CiFQp5S3A9cDrQog6oB54TEqp2QWkK8UppTwghFgJHATqgMellPVaxdmKt4UQMei7Vk4Aj2oaTRNSyjohxBPAt4Aj8JGU8oDGYbWmF7BGCAH6f/efSSk3ahvSz4QQnwPxgE4IkQP8AfgjsFII8SvgFDBTuwj1rhBnvLV+Ptuj7gRWFEWxU6oLSFEUxU6pBKAoimKnVAJQFEWxUyoBKIqi2CmVABRFUeyUSgCK3RJClJthn/WGipC9O/He64QQB5tWmlQUc1LDQBW7JYQol1J6WtM+hRB9gfVSyijTRaUorVNnAIrShBBiqhBit6GAX5IQopfh9QBDTfqfhBBLhRAnhRA6I/ZXLoR4UwiRJoTY1WR/Mw21jtKEENZcjkHpxlQCUJTL7QBGGwr4fQH83vD6H4DvDcXU1gB9jNxfD2CXlHI4+ro7jxhefwW4xfD6NFMFrygdoUpBKMrlQoEVhvouLsBxw+vj0Nd5QUq5UQhxzsj91fBzcbC9wE2Gxz8AiYYSHF+aInBF6Sh1BqAol/s78J6UMhp9TRc3w+utlX02Rq38+UJbPYYfXYbZ4V5CX0k0VQjh3/mQFaVzVAJQlMt5A2cMjx9o8voOYBaAEOJmoEuTfgghBkgpd0spXwGKuLyktKJYhOoCUuyZh6GiY6O/AK8C/xFCnAF2Af0M614DPhdCzAa2oZ/4o6wLx/6zEGIg+jOLzUBaF/alKJ2ihoEqihGEEK5AvaH88xhgiZQyppXt1DBQxWaoMwBFMU4f9LXpHdBf2H3kCttdEEKkApMbZ+AylmHO6A/QdwkpitmpMwBFURQ7pS4CK4qi2CmVABRFUeyUSgCKoih2SiUARVEUO6USgKIoip1SCUBRFMVO/T9HGqxj7mQeJgAAAABJRU5ErkJggg==\n",
|
|
"text/plain": [
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"#####\n",
|
|
"# Plot Correlation\n",
|
|
"#####\n",
|
|
"max_corr_id = np.argmax(corr)\n",
|
|
"min_corr_id = np.argmin(corr)\n",
|
|
"fig, ax = plt.subplots()\n",
|
|
"ax.set_xlabel(\"Lag [ns]\")\n",
|
|
"ax.set_ylabel(\"Correlation\")\n",
|
|
"ax.plot(lag_times/ns, corr)\n",
|
|
"\n",
|
|
"ax.plot()\n",
|
|
"\n",
|
|
"ax.axvline(lag_times[max_corr_id]/ns, ls='--',color='g', label=\"max_corr\")\n",
|
|
"ax.axhline(corr[max_corr_id], ls='--', color='g')\n",
|
|
"ax.axvline(lag_times[min_corr_id]/ns, ls='--', color='r', label='min_corr')\n",
|
|
"ax.axhline(corr[min_corr_id], ls='--', color='r')\n",
|
|
"\n",
|
|
"ax.axvline(calc_lag_time/ns, ls=(0, (5, 5)), color='purple', label='calculated')\n",
|
|
"ax.legend()\n",
|
|
"fig.show();"
|
|
]
|
|
}
|
|
],
|
|
"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.10.2"
|
|
}
|
|
},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 4
|
|
}
|