multilateration.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multilateration in 2D: IoT/LoRaWAN Mass Surveillance\n",
    "\"*Multilateration (MLAT) is a surveillance technique based on the measurement of the difference in distance to two stations at known locations by broadcast signals at known times.*\" - [en.wikipedia.org/wiki/Multilateration](http://en.wikipedia.org/wiki/Multilateration)\n",
    "\n",
    "## Abstract\n",
    "A single motivated individual with several hundred USD and a hobbyist level of competency in electronics and programming would be able to set up a mass surveillance system to track individual LoRa IoT devices on the scale of a small city, provided that the geography of the city is accommodating to the placement of receiving devices (e.g. surrounded by hills).\n",
    "\n",
    "## Locating the source of a signal transmission - intuition\n",
    "\n",
    "When a signal is transmitted by a radio device it propagates through the atmosphere at approximately the speed of light - that is, approximately $v = 3 \\text{e} 8ms^{-1}$. \n",
    "\n",
    "If this signal is received by two receivers (towers), the time difference of arrival ($\\text{TDOA}$) between the first and second tower can be determined. $\\text{TDOA} = t_1 - t_0$ where $t_0$ is the time the signal is received first and $t_1$ is the time the signal is received second. The towers require synchronized high accuracy clocks, which can be achieved in practice with a GPS clock on each tower.\n",
    "\n",
    "The below animation shows a signal propagating from a transmitter \"Tx\" and being received by two towers at different times, allowing the $\\text{TDOA}$ to be calculated.\n",
    "\n",
    "<p align=\"center\"><img src=\"animations/tdoa.gif\"></p>\n",
    "\n",
    "This $\\text{TDOA}$ can be used to determine the *difference* in distance that the transmitter is located from the towers. For example, assume that the transmitter is on a circle of radius $d$ metres from tower 1, where the signal was received first. Then it must also be on a circle of radius $d+v \\times \\text{TDOA}$ metres from tower 0, where the signal was received second. The device must then lie at the intersection of the two circles, giving two possible locations (or one if the circles only touch, or none if the circles do not touch).\n",
    "\n",
    "By iterating over a range of values for $d$ and at each iteration finding the intersection of the two circles, a locus of possible transmitter locations can be determined. This produces a hyperbolic curve on which the transmitter lies. Notably, it is not required to know *when* the signal was first transmitted - as you would with trilateration - and so no communication with the transmitter is required beyond simply identifying the signal.\n",
    "\n",
    "The below animation shows circles of increasing radius around the two towers and the resulting hyperbolic locus of intersections as $d$ is increased. The circle around tower 1 has a radius of $r_1=d$, and the circle around tower 0 has a radius of $r_0=d+\\text{TDOA} \\times v$ (noting $\\text{TDOA} \\times v$ is constant and $\\text{TDOA}$ comes from the previous animation). \n",
    "\n",
    "<p align=\"center\"><img src=\"animations/locus.gif\"></p>\n",
    "\n",
    "With $n$ towers, this process can be repeated between the tower that first received the message and every other tower to produce $n-1$ loci. In general, with three towers the device location can be narrowed down to at least two positions and with four towers to exactly one position. There are some cases, depending on the relative geometry of the towers and transmitter and the error in the timestamping, where this is not the case.\n",
    "\n",
    "<p align=\"center\"><img src=\"animations/loci.gif\"></p>\n",
    "\n",
    "The accuracy of the determined location will depend on the accuracy of the timestamping, the effects of diffraction and obstacles on path length, the relative geometry of the transmitter and towers, etc. The approximation of the problem to two dimensions will also introduce error.\n",
    "\n",
    "## Analysis\n",
    "Multilateration was just explained in an intuitive manner. However, it is convenient to have a set of mathematical expressions that describe the system if the location of the transmitter is to be found. This analysis is inspired heavily by [André Andersen's](http://blog.andersen.im/2012/07/signal-emitter-positioning-using-multilateration) similar work.\n",
    "\n",
    "Consider, in Euclidean $\\mathbb{R}^{2}$ space, a transmitter located at $\\vec{x}$ whose transmission at time $t_0$ propagates at speed $v \\ ms^{-1}$ and is received by a set of $n$ towers. Let the tower that receives this transmission first be at location $\\vec{p}_c$ with receive time $t_c$. Let the remaining $n-1$ towers be located at $\\vec{p}_i$ with transmission receive times $t_i$.\n",
    "\n",
    "For the first tower, $\\vec{p}_c$, we can say that the distance between the transmitter $\\vec{x}$ and the tower is equal to the transmission propagation speed multiplied by the time of flight.\n",
    "$$\n",
    "\\|\\vec{x} - \\vec{p}_c\\| = v(t_c-t_0)\n",
    "$$\n",
    "\n",
    "We can make a similar statement for each of the other towers.\n",
    "$$\n",
    "\\begin{align}\n",
    "\\|\\vec{x} - \\vec{p}_i\\| &= v(t_i-t_0) \\\\\n",
    "&= v(t_i - t_c + t_c - t_0) \\\\\n",
    "&= v(t_i - t_c) + v(t_c - t_0) \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Combining these two expressions we obtain the following, where the only unknown is $\\vec{x}$.\n",
    "$$\n",
    "\\|\\vec{x} - \\vec{p}_c\\| = \\|\\vec{x} - \\vec{p}_i\\| - v(t_i - t_c)\n",
    "$$\n",
    "\n",
    "Expanding this out we obtain a set of $n-1$ expressions, where there are $n$ towers in total. Subscript $x$ and $y$ denote the x and y components of a vector respectively.\n",
    "$$\n",
    "\\begin{align}\n",
    "&\\|\\vec{x} - \\vec{p}_c\\| - \\|\\vec{x} - \\vec{p}_i\\| + v(t_i - t_c) &= 0 &, \\quad i = 1, ..., n-1 \\\\\n",
    "\\Rightarrow \\ & \\sqrt{(\\vec{x}_x - \\vec{p}_{c,x})^2 + (\\vec{x}_y - \\vec{p}_{c,y})^2}\n",
    "              - \\sqrt{(\\vec{x}_x - \\vec{p}_{i,x})^2 + (\\vec{x}_y - \\vec{p}_{i,y})^2} \n",
    "              + v(t_i - t_c) &= 0 &, \\quad i = 1, ..., n-1 \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Each of these equations represents a single hyperbola. These hyperbolas can be plotted individually, producing the same graphs shown above by taking circle intersections, or they can be solved by finding a value of $\\vec{x}$ that minimizes their error. \n",
    "\n",
    "\n",
    "## Implications: IoT and LoRaWAN\n",
    "\n",
    "Multilateration can be used to locate the source of a transmission if \n",
    "1. the transmission can be received in several locations, and \n",
    "2. the receive time can be accurately recorded at each location.\n",
    "\n",
    "Furthermore, if one transmitter's signal contains a persistent and unencrypted transmitter ID then it will be possible to track the location of that device.\n",
    "\n",
    "One emerging IoT technology, LoRa/LoRaWAN, allows these requirements to be satisfied cheaply with hobbyist equipment while also transmitting a pseudo-persistent unencrypted device ID. Every device in a LoRaWAN network is given a 32-bit device address - akin to an IPv4 or IPv6 address - that uniquely identifies it on a network. This address is not globally unique, but is unique within a single LoRaWAN network. This address is assigned by the network operator and can be re-assigned at any time, though typically the same address is used for a substantial length of time, making the address pseudo-persistent - again akin to an IPv4 or IPv6 address. In cases where the device is activated by personalization, the device will always use the same address.\n",
    "\n",
    "Sections 4.3.1.6, 4.3.3, and 4.4.2 of the LoRaWAN [specification](./lorawantm_specification_-v1.1.pdf) specify the encryption of uplink messages (messages sent from an IoT device to the network). Notably, the device address (DevAddr) is never encrypted; in fact it is required that the address be known (unencrypted) in order to decrypt the encrypted portions of the message. Any person listening for LoRa transmissions will receive these messages and will be able to read the device address.\n",
    "\n",
    "LoRa transmissions can typically travel at least 1km, and distances of up to at least 15km can be attained in real-world conditions. As a result, only a low density of receivers need be installed in order to receive all messages over a large geographical area. For example the city of Wellington, New Zealand, is surrounded by hills. A handful of receivers atop, or partway up, some of the hills would be sufficient to cover the entire city. To cover the main portion of the city with three or four overlapping receivers, as required for multilateration, would not be difficult.\n",
    "\n",
    "LoRa receivers with GPS are cheap - the Dragino LoRa/GPS Shield for Arduino comes with a LoRa module and a GPS module with 10ns timing accuracy and typically sells for about 50 USD. This device is additionally easy to configure and use, with an abundance of resources available online.\n",
    "\n",
    "## Summary\n",
    "Given these three points - LoRa messages contain a pseudo-persistent unique ID, LoRa messages can be received over large distances, and LoRa messages can be received and accurately timestamped with cheap hobbyist equipment - there is a very real possibility to create a passive mass surveillance network that tracks individual LoRa devices. This is a severe privacy concern for applications such as vehicle, person (e.g. elderly), and goods tracking.\n",
    "\n",
    "## This notebook\n",
    "\n",
    "This notebook presents a small python script that produces a set of loci of possible transmitter locations given a set of towers with locations and signal receive times. Additionally, the code used to create the animations is given."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "from multilaterate import get_loci\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.optimize import least_squares"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# How many towers. All towers recieve the transmission.\n",
    "num_towers = 4\n",
    "\n",
    "# Metre length of a square containing the transmitting\n",
    "# device, centred around (x, y) = (0, 0). Device will be randomly placed\n",
    "# in this area.\n",
    "tx_square_side = 5e3\n",
    "\n",
    "# Metre length of a square containing the towers,\n",
    "# centred around (x, y) = (0, 0). towers will be randomly placed\n",
    "# in this area.\n",
    "rx_square_side = 25e3\n",
    "\n",
    "# Speed of transmission propogation. Generally equal to speed of \n",
    "# light for radio signals.\n",
    "v = 3e8\n",
    "\n",
    "# Time at which transmission is performed. Really just useful to\n",
    "# make sure the code is using relative times rather than depending on one\n",
    "# of the receive times being zero.\n",
    "t_0 = 2.5\n",
    "\n",
    "# Metre increments to radii of circles when generating locus of\n",
    "# circle intersection.\n",
    "delta_d = int(100)\n",
    "\n",
    "# Max distance a transmission will be from the tower that first\n",
    "# received the transmission. This puts an upper bound on the radii of the\n",
    "# circle, thus limiting the size of the locus to be near the towers.\n",
    "max_d = int(20e3)\n",
    "\n",
    "# Standard deviation of noise added to the\n",
    "# receive times at the towers. Mean is zero.\n",
    "rec_time_noise_stdd = 1e-6\n",
    "\n",
    "# Whether to plot circles that would be\n",
    "# used if performing trilateration. These are circles that are centred\n",
    "# on the towers and touch the transmitter site.\n",
    "plot_trilateration_circles = False\n",
    "\n",
    "# Whether to plot a straight line\n",
    "# between every pair of towers. This is useful for visualising the\n",
    "# hyperbolic loci focal points.\n",
    "plot_lines_between_towers = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "towers:\n",
      " [[ -2006.09826991  -6818.07867897]\n",
      " [   101.31584818   8158.92916584]\n",
      " [-10893.22757653   3893.40549794]\n",
      " [ -8595.14625748   8895.4077718 ]]\n",
      "tx: [ 731.70475905 2184.19298779]\n",
      "distances: [ 9409.38151992  6007.90001384 11749.91315762 11490.45489794]\n",
      "rec_times: [2.50003146 2.50002049 2.50003941 2.50003865]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVcAAAEyCAYAAABOCJPxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xd8VFX+//HXyaQnkEYoSegdQieUhOCKClhBURBWiqJ8UVFXf6xlXRUBkcUKCrsLggvKKiyKiFLEqBhqqFIFQu8lISG9nt8fdxIDJKTNzJ2Ez/PxmMfMnLnlcyG8uTlz7rlKa40QQgjbcjG7ACGEqI4kXIUQwg4kXIUQwg4kXIUQwg4kXIUQwg4kXIUQwg4kXIUQwg4kXIUQwg4kXIUQwg5czS7AXmrVqqUbNWpkdhkCIOkkZCVDnXCzKymTPJ3HkaQjaDRN/ZtiURazSxJOZNu2bZe01sGlLVdtw7VRo0Zs3brV7DIEwLfPwsHVMN75/z7ydT7jYsbBWZjffz7tg9ubXZJwMkqp42VZTroFhChizq45xJ6O5aWIlyRYRaVIuAphtfnsZmbunMndTe5mSMshZpcjqjgJVyGAxMxEXol9hUZ+jXi9x+sopcwuSVRx1bbPVYiy0lrz2vrXSMpK4p+3/xNvN2+zSxLVgJy5Cgdx3nmDF+5fyK+nfuX/df1/tAxsaXY5opqQcBX2Z3GDvByzqyjW/oT9vL/tff5U/08MazXM7HJENSLhKuzP4gF52WZXcZ30nHRe/PVFAjwDmBQ5SfpZhU1Jn6uwP1d3yM0yu4rrvL/tfY5fOc7cfnPx9/Q3uxxRzciZq7A/iwfk50B+vtmVFNpwegOLDixieJvhRNSNMLscUQ1JuAr7c3U3np2kayA5K5nXNrxGE78mPNv5WbPLEdWUdAsI+7N4GM95WeDmaW4twNS4qSRkJDDj1hl4FNQmhI3JmauwP1drgOWaf+YacyKG7458x5j2Y2hbq63Z5YhqTMJV2F9huGaYWkZKdgpTNk2hZUBLnmj/hKm1iOpPugWE/XnUMJ6zUk0tY8b2GVzMuMj0PtNxc3EztRZR/cmZq7C/wnC9YloJv138jUUHFjGs9TDCa1WNeWVF1SbhKuzPw894zkoxZfc5+TlM2DCB2t61eabTM6bUIG4+0i0g7M/kM9cFexcQnxTPjFtn4OPmY0oN4uYjZ67C/jxrGs+Zjg/Xi+kX+feuf3Nr/Vu5tcGtDt+/uHlJuAr7KzxzdXy3wEc7PiInP4fxXcc7fN/i5ibhKuzP3RdcXCEj0aG73Zewj2/iv+HPrf5Mg5oNHLpvISRchf0pBd61IO2Sw3apteadLe/g7+HPmA5jHLZfIQpIuArH8Al2aLjGno5l6/mtPNXxKWq613TYfoUoYJNwVUrNU0pdUErtKdIWqJRao5Q6ZH0OsLYrpdQMpVS8UmqXUqpzkXVGWpc/pJQaWaS9i1Jqt3WdGUom3qx6fIMh7aJDdqW15uMdHxPqG8qgFoMcsk8hrmWrM9f/AP2vaXsZiNFaNwdirO8B7gSaWx9jgH+CEcbAG0B3oBvwRkEgW5cZU2S9a/clnJ2P48L155M/sz9xP2M7jJUrsYRpbBKuWutfgWu/rRgAzLe+ng8MLNK+QBs2Af5KqXpAP2CN1jpRa30ZWAP0t35WU2u9UWutgQVFtiWqioJw1fa9l5bWmn/99i8a1mzIPU3useu+hLgRe/a51tFanwWwPte2tocCJ4ssd8radqP2U8W0X0cpNUYptVUptfXiRcecJYkyqhkCOemQcdmuu9l8bjP7E/fzWPhjuLrINTLCPGZ8oVVcf6muQPv1jVrP1lp31Vp3DQ4OrkSJwub86hvPSSfsupsFexcQ6BnI3U3utut+hCiNPcP1vPVXeqzPF6ztp4D6RZYLA86U0h5WTLuoSvytf7XJJ2+8XCUcSTpC7OlYhrYaKpNgC9PZM1y/BQq+8R8JLCvSPsI6aqAHkGztNlgN9FVKBVi/yOoLrLZ+lqKU6mEdJTCiyLZEVeFnHcSfZL9wXbh/IR4WDwa3HGy3fQhRVjbplFJKfQH8CaillDqF8a3/VGCxUmo0cAJ4yLr4CuAuIB5IBx4F0FonKqUmAVusy03UWhd8SfYkxogEL2Cl9SGqEu9AcPOx25lrek46K46uoG/DvgR6BtplH0KUh03CVWs9tISPbitmWQ08XcJ25gHzimnfClTbSTgTEhK47Tbjj+rcuXNYLBYK+ozj4uJwd3d3SB0ff/wx//73v7FYLNSoUYM5c+bQqlUr22xcKaNrwE59rj+e+JHUnFQeaP6AXbYvRHnJ16lOICgoiJ07dwIwYcIEfH19GT/e/hON5Obm4ur6x4/AiBEjGDduHABff/0148eP57vvvrPdDgMaQ8Jh222viBVHVxDqG0qXOl3ssn0hyksuf3Vy06ZNIzw8nPDwcD766CMApkyZwqxZswB45pln6Nu3LwCrV69m1KhRAKxcuZKePXvSuXNnhgwZQlpaGgBhYWFMmjSJqKgoli5detW+atb84zLRtLQ0bH4hXHALSIiHvFybbjYpM4nNZzbTr1E/29csRAVJuDqxuLg4Fi5cSFxcHBs3bmTWrFns2rWL3r17ExsbC8D27dtJSkoiNzeXdevWER0dzYULF5g6dSoxMTFs376d9u3bM3369MLt+vj4sH79eh566KHr9jljxgyaNm3Kq6++yocffmjbA6rVEvJz4PIxm27255M/k6tz6deon023K0RlSLiaLHn5cg71uY39rdtwqM9tZB44UPhZbGwsgwYNwtvbmxo1ajBw4EDWrVtHREQEW7ZsISkpCV9fXyIiItixYwexsbFER0ezYcMG9u3bR2RkJB07dmThwoUcO3ascLtDhgwpsZ5nn32Ww4cPM3nyZKZMmWLbgw1uaTxfOmjTzW44s4Fgr2BaB7a26XaFqAzpczVR8vLlnH3tdXRmJgC5Z86Qum8vvl5egHEpZ3E8PDwICQlhwYIFREVF0aJFC2JiYjhx4gQtWrRg79699O/fn88++6zY9X18Sr/VybBhw3juueeYO3duBY+uGLWaG88Xf4dWd9lkk/k6n81nN9MrtJd0CQinImeuJrrwwYeFwVpA5+aS+qvxK3/v3r1ZunQpGRkZpKamsmzZMqKjows/e/fdd+nduzfR0dHMnDmTLl2ML3MiIyNZu3YtR44cAYz+00OHDpVaT9Flli9fTsuWLW1ynIU8/cC/IZz9zWabPJB4gMtZl+kZ0tNm2xTCFuTM1US5Z88W2553xbjXVLdu3Rg6dCgREREAPPnkk7Rr1w6A6Oho3nnnHbp3746Xlxdubm6FwVunTh3mzp3LkCFDyM7OBowvwZo3b37Dej788EN++eUX3NzcCAoK4tNPP7XJcV4lpCOc3WmzzW2/sB2AiLoRNtumELagSvrVs6rr2rWr3rp1q9ll3NChPreRe+b6K3ldQ0Jo/lOMCRU5wLoP4McJ8OJR48KCSvpw24fM3zef7Y9sr1S3gLOMNf7555954YUX2L17N0uWLGHgQJkAztkopbZprbuWtpx0C5io9vN/QXl6XtWmPD2p/fxfTKrIAUI6Gc82OntNzEwk0DOw0v2tBWONd+7cydixY3n++ecL39szWHNzrx6W1qhRIxYsWMDgwXIJb1Un4Woiv3vvpd6kibiGhIBSuIaEUG/SRPzuvdfs0uynXgfj+cwOm2wuITOBIM8gm2yrJI4ca9y4cWPatWuHi4v806zqpM/VZH733lu9w/RaXgHGlVqnt9tkc4kZiXadS6DoWOO8vDy6devGLbfcQu/evZk5cyZPPfUU27dvJycnp8Sxxt7e3rz11ltMnz6dv/3tb8AfY41F9SXhKhyvQU84tNq4K0Elf51PzEykiX+TytWzazHETITkU+AXBudbg2934OqxxkDhWOPRo0dfNda4WbNmhWONhw8fftVYY4Ds7Gx69epVuMsbjTUW1YOEq3C8RlHw23+N8a61Kz7wX2tNQmZC5c5cdy2G5c9CTobxPvkkHDgG3p6F+yiOI8Yai6pNOnaE4zWMMp6PravUZtJz08nKy6pcuMZM/CNYC+TnwKE1gOPHGovqQ8JVOF5AI6gZBkfXVmoziRnGdL+VCtfkU8W3ZyYDV4817tGjx3Vjjc+fP0/37t0JDQ0tcaxxhw4diIyM5ODB0i/73bhxI2FhYSxdupTHH3+c9u3bV/zYhKlknKswx/LnYPdX8OIRcK3YUKedF3YyfOVw/nn7P+kV2qv0FYrzQXjxE3j71Yfn91Rsm6Jak3Guwrk17wfZKXBiQ4U3kZCZAFTyzPW218HN6+o2Ny+jXYhKkHAV5mhyC1g84OAPFd5EYqYNugXaD4Z7Z1jvTquM53tnGO1CVIKMFhDmcPeBRr2MIVn9Kza1oU36XMEIUglTYWNy5irM06KfcWeCS/EVWt3iYgEgN9+2dzYQwhYkXIV5Wt0DysUY81oBBZe9FnQPCOFMJFyFefxCoXlf2PE55OWUe/UgLwlX4bwkXIW5uoyC1PNwcFW5Vy04c72UccnGRQlReRKuwlzN7oAaIbDtP+VeNcQ3BIDTqadtXJQQlSfhKsxlcYXOwyE+Bi4fL9eqAZ4B+Hn4cTT5qJ2KE6LiJFyF+ToNN553FD/JyY00qtmII8lHbFyQEJUn4SrM518fmt9h/WKrfMOqWgW2Yn/CfvLy8+xUnBAVI+EqnEOXUZByFg6V74qtDsEdSM9N53DyYfvUJUQFSbgK59C8H/jWLfcXWx2CjdvG7L642w5FCVFxEq7CORR+sbUGkoqZpaoE9WvUx9/Dn12XdtmxOCHKT8JVOI9Ow41bv+z4vMyrKKXoENyBredkeknhXCRchfMIaAjNbjNGDZTji62eIT05kXKCk1fKfsYrhL1JuArn0uVRuHIa9i4tfVmr6FBj9v91Zyp32xghbEnCVTiXlndB7bawdmqZz14b1GxA/Rr1WXdawlU4DwlX4VxcXODWV4ypCHcvLvNqt4Tdwuazm0nLSbNjcUKUnYSrcD6t7oG67WHtP8o8W1bfRn3Jysvi55M/27k4IcpGwlU4H6Xg1lfh8jHYWba5XjsEd6COdx1WHS3/7FpC2IOEq3BOLfpBaFf4ZSpkl/6rvotyoV+jfqw/s57krGQHFCjEjUm4CuekFPR7C1LOQOz7ZVqlf6P+5Obn8tOJn+xcnBClk3AVzqtBD2g3GDZ8BImlTysYXiucBjUasOzwMgcUJ8SNSbgK53bHm+DiCj/8vdRFlVLc3/x+tp3fxrHkY/avTYgbkHAVzq1mCPQeD79/Z0yoXYoBTQdgURa+jv/aAcUJUbJqG66X07JJSs82uwxhCz2fhoDGsOrlUodmBXsHc0vYLSyLX0ZOfvlveiiErVTbcD2VlEHXyT8ycl4ci7eeJDld/qFVWa4e0P9tuHQQNn5c6uKDWgwiMTORmBOln+kKYS/VNlybBfsyOroxhy+m8uKSXXSZvIYR8+L4Iu4ECalZZpcnyqtFf+Pigl+mwqX4Gy7aK7QX9WvUZ+G+hQ4qTojr2T1clVLHlFK7lVI7lVJbrW2BSqk1SqlD1ucAa7tSSs1QSsUrpXYppToX2c5I6/KHlFIjS9uvl7uFV+5sTeyLt/LtuChG92rMsUtpvPL1biLe+pGHZ29k/oZjnEvOtN/BC9tRCu5+zziL/XYc5OeXuKiLcmFYq2HsvLiTPZf2OLBIIf6gtNb23YFSx4CuWutLRdqmAYla66lKqZeBAK31S0qpu4BngLuA7sB0rXV3pVQgsBXoCmhgG9BFa325pP127dpVb9169RyfWmv2nb3Cqj3nWLnnHPEXUgHo3MCfO8Pr0T+8LvUDvW149MLmdiyEZU/BXe9CtydKXCw1O5Xbl9zOn+r/ianRUx1YoKjulFLbtNZdS13OpHA9APxJa31WKVUP+EVr3VIp9W/r6y+KLlfw0Fr/n7X9quWKU1y4Xiv+Qkph0O49cwWAtiE1uTO8Lv3D69Gstm+Fj1vYidbw+SA4sQme2mjMAVuCf8T9gy9//5LVD66mtndtBxYpqrOyhqsj+lw18INSaptSaoy1rY7W+iyA9bngJz8UKDrj8SlrW0ntldKsdg3G9WnO989G8+tfb+Vvd7XC3dWFd384yO3vr+WO99fy/g8H2HfmCvb+T0iUkVJw74egXGDp/91wWsJhrYaRp/P44vcS/w8Wwm5cHbCPKK31GaVUbWCNUur3GyyrimnTN2i/emUjvMcANGjQoFxFNgjyZkzvpozp3ZSzyRmstp7RfvxzPDN+iqdhkDf9w+tyZ3g9OoT5oVRxJQmH8G9g9L8uHQPr3odbXix2sfo163N7w9tZ9PsiHgt/jBruNRxcqLiZ2f3MVWt9xvp8AVgKdAPOW7sDsD5fsC5+CqhfZPUw4MwN2q/d12ytdVetddfg4OAK11zPz4tRUY1Z9H89iXv1dt5+oB0Ng3yYG3uUgTPXEzn1JyZ8u5fNRxLIy5czWlN0GGJcGvvLVDgZV+Jio9uNJiUnhUUHFjmwOCHs3OeqlPIBXLTWKdbXa4CJwG1AQpEvtAK11i8qpe4GxvHHF1oztNbdrF9obQMKRg9sx/hCK7GkfZelz7W8ktNz+HH/eVbuOcevhy6SnZtPLV8Pbm9dm75t6xDZtBaebhab7lPcQGYy/KsXoGDsOvCsWexiY38cy/6E/awatAovVy/H1iiqHaf4Qksp1QTjbBWMLoj/aq3fUkoFAYuBBsAJ4CGtdaIyftf+GOgPpAOPaq0Lhm89BvzNuq23tNaf3mjf9gjXolKzcvn59wus2nuOtQcukpqVi7e7hVtaBNO3bR36tKyDn7eb3fYvrE7Gwbz+0O5BeGB2sYtsO7+NUatG8XK3l/lz6z87uEBR3ThFuJrJVuGakJDAbbfdBsC5c+ewWCwUdDnExcXh7u5OVm4eGw8nsGbfedbsO8+FlCwsLorujQPp26YOd7StS6h/5c6Ypk2bxqeffoq7uzu1a9dm3rx51K9fv/QVbwa//AN+mQIDZkKnR4pdZOTKkZxOPc3KB1biZpH/9ETFSbja4cx1woQJ+Pr6Mn78+BKXyc/X7DqdzA97z7Fm33kOWcfStg2pyR1t6tC3TV1a16tR6hdiubm5uLr+8X3jTz/9RM+ePfHy8uKjjz5i06ZNLFwoVyABkJ8Hnw00zmJHr4F67a9bZN3pdTz545NM6DmBQS0GmVCkqC6caShWtTVt2jTCw8MJDw/no48+AmDq1LfZsPy/vNi/Fa2OLqH+pvd55c5WJB/ayt+ff4q7ZsQSPnoq9Vt1pGXb9gwePJi0NGOm/bCwMCZNmkRUVBRLl159a+k+ffrg5WWc/fbo0YNTp0459mCdmYsFBs0Dr0BYPBwykq5bJCokivCgcGbvmk1OGe/LJURlSLhWUFxcHAsXLiQuLo6NGzcya9Ysdu3aRe/evYmNjQVg+/btZKWnMDqqIbf6JzJ57IP87dZ6JG74Hx4D3iDr3reJTfCm3+Mvs2rPWTTg4+PD+vXreeihh0rc99y5c7nzzjsddKRVhG8wDJ4Pyadg6djrLo9VSjGu0zjOpJ3h60MyHaGwPwnXCoqNjWXQoEF4e3tTo0YNBg4cyLp164iIiGDLli0kJSXh6+tLREQEO3bsIDY2lrvu6EPtjBPkJpzEe9UELN+8SPaBX9l/8DBjP9/O+SuZrM9vzqfrj3I8ofj7Rs2fP5/du3fzwgsvOPiIq4D63aDfFDi40hj/eo3IkEg61+7M7F2zycyVOSWEfTniIoIq6fsj3zN9+3TOpZ2jrk9dnuv83FWfl9RX7eHhQUhICAsWLCAqKooWLVoQExPDiRMnaNGiBXv37qV///589tlnhevk5uWz5dhl+n/iyqUsxZvL9/Hm8n00q+1Ln1a16dOqNl0aBhCz5gemTZvG2rVrcXd3t+vxV1ndxhh9rz9Nhjrh0LJ/4UcFZ6+PrX6MRQcWMbJtqfP/CFFhcuZajO+PfM+EDRM4m3YWjeZs2lkmbJjAwcsHC5fp3bs3S5cuJSMjg9TUVJYtW0Z0dHThZ++++y69e/cmOjqamTNn0qVLFwAiIyNZu3YtR44cASAtLY2jRw7Ts2kQNTxdWfZ0L9b+9U+8cW8b6vl58un6ozw8exOtn5zJkJFPMO7t2Vi8ih/PKTAuj73vI6jXAb4aDRf2X/VxRN0Ietbrydzdc0nPSTepSHEzkHAtxvTt08nMu/rXxsy8TDaf3Vz4vlu3bgwdOpSIiAh69OjBk08+Sbt27QCIjo7m/PnzdO/endDQUNzc3AqDt06dOsydO5chQ4bQoUMHIiMjOXjw4FX7ahjkw6NRjflsdHd2vN6Xfz3Smaz180lPT+O5J4ZTt0krQjvewqxf4jlwLkXmPbiWuzc8/F9w94EvHoa0hKs+HtdpHJezLrNwv4y2EPYjQ7GK0X5+e/T1UxegUOwauauypVVYfr5m9+lkYn6/wE+/n2fPaWMmr1B/L3q3COaWFrWIbFaLmp4yjhOAU1vh07uMvtjhS6HI+NZnYp5h24VtrHxgJX4efiYWKaoaGedaiXDtu6QvZ9POXtdez6cePzz4Q2VLs5nzVzL5+fcL/PT7BTYcTiA1KxeLi6JzA39uaRFM7xbBhIf44eJyE08y89siY4KXLo/CPR8Y3QbAgcQDPLT8IUa1HcULXeXLQVF2Eq6VCNeCPteiXQOeFk8mRE7g7iZ326pEm8rJy2fHiSTWHrzArwcvsft0MgCBPu5EN69F7+bBRLeoRe0aniZXaoI1b8D6D42RBD2fLmx+dd2rrDq6iuX3LyfEN8TEAkVVIuFaySu0ihst4KzBWpxLqVmsO3SJtQcvEnvoIpdSjTvhtqlX09qFEEyXhgG4u94E3e75efC/UbB/uTEWts0AAM6lneOepfdwR8M7eDv6bXNrFFWGhKudJ26pSvLzjdvbrD14kbUHL7L9+GVy8zVebha6NgogqlktoprWok1ITSzVtQshJwMWDIAzO2Hkt9CgBwAfbvuQuXvmsuieRbQJamNykaIqkHCVcC1RSmYOGw8nsOFwAuvjLxXOf+Dn5UaPJoFENatFZNMgmgb7Vq9JwdMSYO4dkJFozEFQqzkp2Snc9fVdtAxoyZy+c6rX8Qq7kHCVcC2zCymZRtjGJ7D+8CVOXc4AoHYNDyKbBhFpDduwgGpw88bEI/DJHcYwrcd/BN/aLNy/kKlxU5l12yyiw6LNrlA4OQlXCdcKO5GQzobDl1h/OIGNhy8V9teG+nsR0SiAiMaBdGsUSLPaVfTM9tQ2+M/dULsVjPqeHIs7A5YNwMPiwZJ7l2BxkQnPRckkXCVcbUJrzcHzqWw8fIktxy4TdyyRiylZAAR4u9G1kRG0EY0DaRtSEzdLFfmC7MBK+HIYNL0NHv4vq0/9zPi142VKQlEqCVcJV7vQWnM8IZ24o4nEHUtky7FEjicYl5F6u1vo1MCfiEaBRDQKpH2YHzWc+YKGbfNh+bPQ9gH0A3MY+cNjHL9ynO/v/x5fd7mtuiheWcNVJm4R5aKUolEtHxrV8mFwhHEnhAtXMo2gPZpI3LHLTI85hNbGeP3mtX3pVD+Ajg386dTAn+a1azjPiIQuIyEzCda8jvLy56XuL/LwiqHM3jVbLiwQlSZnrsLmkjNy2HHiMjtPJhU+ktKNCap93C20C/OjU4MAOoT50y7MjxA/T3P7bn+cAOs+gOjxvOaZzXdHvuObAd/QsGbDwkXKcrsfcXOQbgEJV6ehteZYQjo7T15mxwkjbPeduUKu9bbkAd5uhIf60TbEj/DQmoSH+NEg0Ntxl+1qDd/9Bbb9h0t9/sbdJ7+mW71ufNTno2IXL8vtfkT1Jd0CwmkopWhcy4fGtXy4v1MYAJk5eew7e4W9Z66w93Qye84kM3fdEXLyjMCt4eFK65CatK5bg1b1atKybg1a1qmBj4cdfmSVgrvfh8xkav00hTE9h/PhyV/YcGYDkSGRN1x1y5YtjB49mri4OPLy8ujWrRuLFi0iPDzc9nWKKkXCVZjC081C5wYBdG4QUNiWnZvPwfMp7D2TzJ7TV9hzJpkl206Rlp1XuEyDQG9a1KlBoyBvGgZ50yDIh4aB3oQGeFVupIKLBe6fDZnJDN+0kK9atmda3DSW3LcEV5eS/5lERERw33338fe//52MjAweeeQRCVYBSLeAcHL5+ZrTSRnsP3uF38+lcOBcCgfPp3AiMZ2s3D/uk2VxUYT4e9Iw0IcGQd40CvKmQaCPEcCB3mU/481Og88eICZxD3+pHcjoOk8T8Gs7UhOz8A30oOeApvx35b+u6hbIzs4mIiICT09PNmzYgMUi42SrM+kWENWCi4uifqA39QO96du2bmF7fr7mQkoWxxPSOJ6YzomEdOtzGit2ny38Aq1ALV8PGgZ50zDQ2xq+Rgg3DPQm0Mf9jy/U3H3gz4vps2AA3TPO8N/TcxhyZQKe+JCamMXPC38nIT0V31Z/DNVKTEwkNTWVnJwcMjMz8fHxccifjXBuEq6iSnJxUdT186SunyfdmwRd93lyRo41cNM4npBe+HrjkQSW7jxN0V/YfD1caRBY0M3gTcNAH0K6zWb0d4MZG5bFvoaf0fnwWABys/M5eSCRhq3+CPoxY8YwadIkjh49yksvvcTHH39s9+MX9qG1Jis3n4zsPNKyc0nPziM1K5c06yM1K6/0jVhJuIpqyc/LjXZhfrQLu/4uA5k5eZy6nM7xBONxIjGd4wlpHDifQsz+C2TnGd0N/rzMHVcmsyZ4H0lpP3I0rRfeWnE2N5v0vefw/OkQe9d+x6W0XBpE9KVNlAuPPtCPb1f8wD3977i5Jym3s7x8TWZOHpk5eaRnG4+07FwjFLNyycjJIy0rj3RrQP7xWR4ZObnGc5EATc/OJT0rj/ScPPLybdNVKn2uQhSRl685fyWTiylZLJq1k7SsE2wLn0ajnBxqxY9lq25Fhitke7qQkplb4naUgpqebvh5GY+aXq54u7vi7W7B290VH3eL8drDeO3lXvBswcPVgrurws3igpvFBXdXF9yLvHazGJ+5W1wcFuBaa/I15GtNbp4mOy+fnIJH7jXv8/LJztVXv8/T5OQAhMGnAAAfkUlEQVQWvM63BmOR51wjKLMK2nKLfm68zirSllvOAHR3dbH+mRf8HRR57eGKt5sFb4+r233cXfFyt+Dr6Yqvhys+7sZzw1o+0ucqRHkZX4x5EeLvhdfAVvy8EFyPDeCHpt8yofaH9D7zKuEP3UeL7nXJzMnjcno2yRk5JKfnGM8lPK5k5JCQml14llVwRlVZri4KV4vCRRkPhRHsSilcrnlWYCxjzeN8rcnLLwhOTV6+RmvIs77Pz7cuozX2OgdzdVF4ulnwdHPBw9V4Nt4br/293PB0s+BR0H7VMsazl5sFHw8jCH2uDU8PC95uFlxNmPNCwlWIErTobvSreiyzsDtlG9MDNcuyJxNQPwKoi6ebhXp+XtTz86rQ9vPzNZm5edf9ipqda5zdFT3Ty87NJyfvj7PBrNyiZ4r6qjNLrY3A1PzxPr+gzboMGP+RFISvxeWPgC54r5TC4kKRduMzN1frWbTlj7NrN1fjvauLi/VzVXi2bZxxX30mbgSliymh5yjSLSBEGexL2MfQ74bycGY+r1xONu4mG9rF7LKECco6FKv6/rchhA21CWrD4JaD+dLThd99/GDBQDi5xeyyhBOTcBWijMZ1Goe/pz9TGrUh3ycIPrsfTmwyuyzhpCRchSgjPw8//tL5L+xI3MfyW5+HGnXgswfg2HqzSxNOSMJViHIY0GwAHYI78P6+eVwZtgj8QmHhg3D0V7NLE05GwlWIcnBRLrza/VWSspL4OH4xjPoe/BvCwsFw+GezyxNORMJViHJqHdSawS0Gs+jAIvZnJcCo7yCwCXzxMMT/aHZ5wklIuApRAeM6jcPfw5/JmyaT5xUAI5dDrebwxVA4uNrs8oQTkHAVogL8PPwY33U8uy7t4qtDX4FPEIz4Fmq3hi//DPuWmV2iMJmEqxAVdE+Te+hetzsfbvuQSxmXwDvQCNjQzvC/UbDzC7NLFCaScBWigpRSvNrjVTLzMpm2ZZrR6OUPj3wNjaLhm7Gw5RNzixSmkXAVohIa+zXm8XaPs/LoSjac3mA0evjCsMXQ4k74/v/B+unmFilMIeEqRCWNbjeahjUbMnnzZDJzM41GN08Y8hmED4I1r8NPk7Hb1FLCKUm4ClFJHhYPXuvxGidTTjJn95w/PrC4wQNzoNNw+PUdWP03CdibiISrEDbQvV537mlyD/P2zONI0pE/PnCxwH0fQY+nYNMsWP4s5Fd+Hlfh/CRchbCR8V3H4+3qzaRNk7hqKk+loN8U6P0ibF8AXz8BeTklb0hUCxKuQthIkFcQz3d5nq3nt7Ls8DXjXJWCPq/CHRNhz1ewaDjkZJpTqHCIKhOuSqn+SqkDSql4pdTLZtcjRHEeaP4AnWp34r2t73E58/L1C0Q9B3e/BwdXweeDIDPZ8UUKh6gS4aqUsgAzgTuBNsBQpVQbc6sS4nouyoXXerxGanYq7297v/iFIh6HQZ/AyU3wn3sg9aJjixQOUSXCFegGxGutj2its4EvgQEm1yREsZoHNGdE2xF8E/8NW8+VcKuhdg/C0C/h0iGY1w+STji2SGF3VSVcQ4GTRd6fsrZdRSk1Rim1VSm19eJFORsQ5hnbYSyhvqFM3DSR7Lzs4hdqfgeMWAbpl2BuP7iw37FFCruqKuFa3M3ZrxswqLWerbXuqrXuGhwc7ICyhCiel6sXr3Z/laPJR5m7Z27JCzboDo+uBJ0P8/rLfbmqkaoSrqeA+kXehwFnTKpFiDKJDovmzsZ3MmfXnKvHvl6rTlsYvRq8AmDBfRAf47gihd1UlXDdAjRXSjVWSrkDDwPfmlyTEKV6KeIlvN28mbBxAvk6v+QFAxrBY6shsCn8dwjs+dphNQr7qBLhqrXOBcYBq4H9wGKt9V5zqxKidEFeQYzvOp4dF3aw5OCSGy9co45xV4OwCFjymMyoVcVViXAF0Fqv0Fq30Fo31Vq/ZXY9QpTVgKYD6F6vOx9s+4DzaedvvLCXPwz/Glr0N2bUWvuOzEdQRVWZcBWiqlJK8XqP18nJz+HtuLdLX8HNy5hRq/3D8PNkWPUK5N+gS0E4JQlXIRygQc0GPNnhSWJOxPDj8TLcxNDiBgP/CT2ehs3/NOYjyM2yf6HCZiRchXCQEW1H0CqwFVM2T+FK9pXSV3BxgX5vwe1vwp4lsPAhyCzDesIpSLgK4SBuLm5MiJxAQmYCH277sGwrKQW9/gID/wXH18N/7oaUUvpthVOQcBXCgdoGteWR1o/wv4P/Y9v5bWVfseNQGLoIEg7D3DuMZ+HUJFyFcLCnOz5NqG8ob258k6y8cvSjNr8dRi2H7FQjYE+XI5yFw0m4CuFg3m7evNbjNY4mH2XOrjmlr1BUaBcYvQbcfY0ZtQ6V4csxYQoJVyFMEBUaxd1N7mbunrnEX44v38pBTY2ADWoGXwyBnV/Yp0hRKRKuQpjkxYgX8XXzLf3S2OLUqAOjvoeGUfDNWFj3gVxs4GQkXIUwSaBnIC9GvMhvF39j0YFF5d+AZ0348xIIfxB+nCAXGzgZCVchTHRPk3voWa8n07dP51zaufJvwNXduH13z3HGxQZfPSYXGzgJCVchTKSU4vWer5OXn8fkTZOvvmtsWRVcbNB3MuxdKvfmchISrkKYLKxGGOM6jWPtqbWsPLqy4huKfMY4iz2xET69G1IqcCYsbEbCVQgn8EjrR2hXqx1T46aSmJlY8Q21HwzDFsPlo/DJHXDxgO2KFOUi4SqEE7C4WJgYOZGUnBSmxk2t3Maa3WaMJMjLMi42OLbeNkWKcpFwFcJJNAtoxpj2Y1h5dCW/nPylchsL6WiMhfWtC58NhN2lTNQtbE7CVQgn8nj44zQPaM6kjZNIyU6p3MYCGhr35gqLgK9Gw/rpMhbWgSRchXAibhY3JkZO5FLmJd7b+l7lN+gVAI98DW3vhzWvw4q/Qn5e5bcrSiXhKoSTCa8Vzsg2I/nq0FdsPru58ht084RB84zRBFvmwKLhkJ1e+e2KG5JwFcIJPdnxSRrUaMCEDRNIz7FBELq4GONg73wHDqyA+fdC2qXKb1eUSMJVCCfk5erFhMgJnEo9xcc7P7bdhruPMe7PdX6PzAtrZxKuQjipiLoRDG4xmM/3fc5vF3+z3YZb3wsjl0NGkhGwJ7fYbtuikISrEE7s+S7PU9u7Nm+sf4PsvGzbbbh+N3j8R/CoaXQR/P697bYtAAlXIZyar7svr/d8ncPJh5mzu5wTa5emYF7YOm3gyz9DnI23f5OTcBXCyfUO6809Te7hk12fcCDRxpez+gbDyO+g5Z2wYjz88JpMW2gjEq5CVAEvRbxETY+avLHhDXLzc227cXdvGPI5RDwOG2YYFxzkZNp2HzchCVchqgB/T39e6f4KexP28tm+z2y/AxcL3PUu3P4m7P0aPn8A0isxgYyQcBWiqujXsB996vdh5s6ZHEs+ZvsdKAW9/gKD5sKpLTCvP1w+bvv93CQkXIWoIpRSvNrjVdxd3Ct2362yavcgDF8Kqefgk9vlFt4VJOEqRBVS27s2f434K9vOb+N/B/5nvx016mWMJHDzNCbelqFa5SbhKkQVM7DZQHrU68H7297nbOpZ++0ouCU8HvPHUK1N/7TfvqohCVchqhilFG/0fAONZuKmiRW771ZZ+dY2hmq1uhtWvQwrXpRZtcpIwlWIKiisRhjPdX6OdafXsfzIcvvuzN0bBi8w7jAb92/jLDY7zb77rAYkXIWoooa2Gkqn2p2YGjeVi+kX7bszF4txh9m73oVDq+HTu+QGiKWQcBWiinJRLkyMnEh2Xrb9uwcKdHsCHv4CLh0yRhKc32f/fVZREq5CVGGN/BoxruM4fjn5CyuOrnDMTlv2h0dXQF4OzOsHh39yzH6rGAlXIaq44W2G0z64PW/Hvc2lDAdNgB3SEZ6IAb/6sPAh2L7AMfutQiRchajiLC4WJkVOIiMngymbpzhux35h8NgqaHwLfPsMxEyUSV+KkHAVohpo4t+EJzs+yZrja1h9bLXjduxZE4Ytgs4jIfY9+PpxmfTFSsJViGpiVNtRtA1qy1ub3iIx04GTrljc4N7pcPsE2PMVfDZQJn1BwlWIasPVxZVJUZNIyUnh7c1vO3bnSkGv5+HBT+H0dmMkwU1+fy4JVyGqkeYBzRnbfiyrjq0i5niM4wsIf8B6f67LRsCe2OT4GpyEhKsQ1cxj7R6jdWBrJm2aRFJmkuMLaNDduD+XVwDMv8/oKrgJSbgKUc24ubgxKWoSyVnJTN0y1ZwigpoaARvaGZY8BrHvgyMucnAiEq5CVEMtA1vyRPsn+P7I9/xy8hdzivAOhOHfQPiDEPMmLH/WuPDgJmG3cFVKTVBKnVZK7bQ+7iry2StKqXil1AGlVL8i7f2tbfFKqZeLtDdWSm1WSh1SSi1SSrnbq24hqosn2j1Bi4AWTNw4keSsZHOKcPOEB+ZA9HjjQoOFD0GmSbU4mL3PXD/QWne0PlYAKKXaAA8DbYH+wCyllEUpZQFmAncCbYCh1mUB/mHdVnPgMjDaznULUeW5WYzugcTMRN7Z8o55hbi4wG2vwYCZcCwW5va9KW4fY0a3wADgS611ltb6KBAPdLM+4rXWR7TW2cCXwACllAL6AEus688HBppQtxBVTpugNjwW/hjLDi8j9lSsucV0esS4fUzKWfjkNji5xdx67Mze4TpOKbVLKTVPKRVgbQsFThZZ5pS1raT2ICBJa517Tft1lFJjlFJblVJbL1608xRsQlQRYzuMpZl/MyZsnEBKdoq5xTTubdzdwN0X/nN3tR5JUKlwVUr9qJTaU8xjAPBPoCnQETgLvFewWjGb0hVov75R69la665a667BwcHlPh4hqiN3izuToiZxKeMS7219r/QV7K1WcyNgC0YS/PpOtRxJ4FqZlbXWt5dlOaXUHOA769tTQP0iH4cBZ6yvi2u/BPgrpVytZ69FlxdClEF4rXBGtR3FvD3z6NuwL5GhkeYW5BMEI5YZE778NBkSjsC9H4Krh7l12ZA9RwvUK/L2fmCP9fW3wMNKKQ+lVGOgORAHbAGaW0cGuGN86fWtNmYA/hl40Lr+SGCZveoWorp6quNTNPZrzISNE0jLcYLbtLh6wP3/hj/9DX77L3x2f7Wak8Cefa7TlFK7lVK7gFuB5wG01nuBxcA+YBXwtNY6z3pWOg5YDewHFluXBXgJeEEpFY/RBzvXjnULUS15WDyYGDmRc2nneH/r+2aXY1AK/vQSPPAJnNpSreYkUA65NYQJunbtqrdu3Wp2GUI4nXe2vMOCfQv4pO8ndK/X3exy/nBiE3w5DHQ+DPkcGvUyu6JiKaW2aa27lracXKElxE1mXKdxNKzZkDc2vEF6TrrZ5fyhQQ/jklmfYFgwEHZ+YXZFlSLhKsRNxsvVi4mREzmTeoYPt39odjlXC2wCo3+Ahj3hm7EQM6nK3t1AwlWIm1DnOp0Z1noYX/z+BVvPOVn3mVcAPPI1dBoOse/CV6MhJ8PsqspNwlWIm9SznZ4lzDeM1ze8Tkauk4WXxQ3u+wjumAh7l8L8eyH1gtlVlYuEqxA3KW83byZGTeRkykk+2vGR2eVcTymIeg4GL4Bze4xLZi/sN7uqMpNwFeImFlE3giEth/D5vs/Zfn672eUUr8198Oj3kJtlTPoSb8IdFipAwlWIm9wLXV4gxDeE19a/5nzdAwVCuxiXzPrVN6Yt3DrP7IpKJeEqqoyEhAQ6duxIx44dqVu3LqGhoYXvs7OzHVZHZmYmDz74IM2aNaNnz56cOHHCYfu2B283byZFTeJEyglmbJ9hdjkl868Po1dD0z7w3fOw+lXIzzO7qhJJuIoqIygoiJ07d7Jz507Gjh3L888/X/je3d1+86fn5uZe9X727NnUrVuX+Ph4nn76aV555RW77dtRIupG8HDLh1m4fyHbzm8zu5ySedSAoV9CtzGw8WNYNByyneBS3mJIuIpqYdq0aYSHhxMeHs5HHxlfzkyZMoVZs2YB8Mwzz9C3b18AVq9ezahRowBYuXIlPXv2pHPnzgwZMoS0NOMfalhYGJMmTSIqKoqlS5deta9ly5YxcuRIAAYPHszq1asdcYh293yX5wu7B5zq4oJrWVzhrnfgzmlwcCXM6w9XnG8uJwlXUeXFxcWxcOFC4uLi2LhxI7NmzWLXrl307t2b2Fhjgujt27eTlJREbm4u69atIzo6mgsXLjB16lRiYmLYvn077du3Z/r06YXb9fHxYf369Tz00ENX7e/06dPUr29M4Obu7o6Pjw9JSSbcZdXGCroHTqacZMYOJ+4eKND9/4yz2MQjMOc2OPub2RVdRcJVOL1vdpwmaupPNH75e6Km/sQ3O05f9XlsbCyDBg3C29ubGjVqMHDgQNatW0dERARbtmwhKSkJX19fIiIi2LFjB7GxsURHR7Nhwwb27dtHZGQkHTt2ZOHChRw7dqxwu0OGDCm2nuLm4zBumFH1RdSNYFirYSzcv5At56rAnQJa9IPHVhnDtubdCQdWml1RIQlX4dS+2XGaV77ezemkDDRwOimDV77eze9nrxQuU9LkQx4eHoSEhLBgwQKioqKIjo4mJiaGEydO0KJFC7TW9O/fv7Dfdt++fcyePbtwfR8fn2K3GxYWxsmTxk0zsrOzSUtLw8/Pz3YHbbLnOj9H/Rr1eX39687dPVCgbjt44idjEu4vhsLGWU4x+baEq3Bq76w+QEbO1d8IZ+Tksf5wQuH73r17s3TpUjIyMkhNTWXZsmVER0cXfvbuu+/Su3dvoqOjmTlzJl26dAEgMjKStWvXcuTIEQDS0tI4dOhQqTXdd999zJ8/H4DFixcX9uVWFwXdA6dTTzvf3AMlqVEXHl0Bre6G1a/A9y+YfhtvCVfh1M4kFT/uMiXzj3843bp1Y+jQoURERNCjRw+efPJJ2rVrB0B0dDTnz5+ne/fuhIaG4ubmVhi8derUYe7cuQwZMoQOHToQGRnJwYMHS61pzJgxnD17lmbNmvHxxx8zZcoUGxypc+lSpwt/bv1nvvj9i6rRPQDg7gODPzOu6to6zxgPm2FeX7jM5yqcWtTUnzhdTMCG+nux/uU+JlR088jIzeDBbx8kT+fx9X1f4+3mbXZJZbfjc1j+FwhsDMMWGbNt2YjM5yqqhb/2a4mXm+WqNi83C3/t19Kkim4eXq5eTIqaxJnUM7y/zUnuXFBWnR6BEd9A2kVjJMHxDQ4vQcJVOLWBnUJ5+4F2hPp7oTDOWN9+oB0DOxV7d3VhY53rdOaRNo+w6MAiNp/dbHY55dOol3HJrHcgzL/P4ZNvS7eAEOKGMnIzeGj5Q+Tm5/LVfV/h41b8KAqnlXEZFo+Ao79C9P+DW/8OLhU/r5RuASGETVzVPeAsNzYsj4LJtzuPhNj34H8jIdv+Q8wkXIUQpepUuxMj2oxg8cHFbDq7yexyys/iBvdOh75vwf7l8J+7IOWcXXcp4SqEKJNxncbRqGYjXl//OqnZqWaXU35KQeQ4GPoFXDwIc/rA2V12252EqxCiTDxdPZkUNYnz6ed5b9t7ZpdTcS3vNKYuBGPSl99X2GU3Eq5CiDLrWLsjI9qMYMnBJWw44/jhTTZTcMlscEv4chisn2HzS2YlXIUQ5fJ0x6dp7NeYNza8UTW7BwrUqAujvjduI7PmNVj+LOTabtJ1CVchRLl4unoyOWoyF9Iv8O7Wd80up3LcveHB/0Dvv8L2BfD5A5CeaJNNS7gKIcqtfXB7RrYdyVeHvmL96fVml1M5Li7Q5+9w/7/h5GaIedMmm5WLCIQQFZKVl8Xg5YNJy0lj6YCl1HCvYXZJlXdyCwS3AM+Sp5CUiwiEEHblYfFgctRkLmZc5J0t75hdjm3Uj7hhsJaHhKsQosLaBbfj0baPsjR+KbGnYs0ux6lIuAohKuWpjk/RzL8ZEzZO4Er2ldJXuElIuAohKsXd4s7kqMkkZCQwLW6a2eU4DQlXIUSlta3VlsfCH2PZ4WX8eupXs8txChKuQgibGNthLM38m/HmhjdJzko2uxzTSbgKIWzC3eLO5F6TSchMYNoW6R6QcBVC2EzboLaMbjeabw9/yy8nfzG7HFNJuAohbGps+7E0D2jOxI0Tb+ruAQlXIYRNuVncmBw1mcuZl/lH3D/MLsc0Eq5CCJtrE9SGx9s/zvIjy/n5xM9ml2MKCVchhF2MaTeGlgEtmbjp5uwekHAVQtiFm8WNyb0mk5SZxJTNU8wux+EkXIUQdtMqsBVj2o9hxdEVxByPMbsch5JwFULY1ePtH6d1YGsmbppIYqZtJqKuCiRchRB25eZidA9cyb5yU3UPSLgKIeyuRUALnurwFKuPrWbVsVVml+MQlQpXpdRDSqm9Sql8pVTXaz57RSkVr5Q6oJTqV6S9v7UtXin1cpH2xkqpzUqpQ0qpRUopd2u7h/V9vPXzRpWpWQhhjkfDHyU8KJy3Nr3FpYxLZpdjd5U9c90DPABcNQ2OUqoN8DDQFugPzFJKWZRSFmAmcCfQBhhqXRbgH8AHWuvmwGVgtLV9NHBZa90M+MC6nBCiinF1cWVyr8mk56QzedNkqustpgpUKly11vu11geK+WgA8KXWOktrfRSIB7pZH/Fa6yNa62zgS2CAUkoBfYAl1vXnAwOLbGu+9fUS4Dbr8kKIKqapf1PGdRpHzIkYVhxdYXY5dmWvPtdQ4GSR96esbSW1BwFJWuvca9qv2pb182Tr8kKIKmhEmxF0q9uNtJw0s0uxK9fSFlBK/QjULeajV7XWy0parZg2TfFhrm+w/I22df1OlRoDjAFo0KBBCaUJIcxkcbHwSd9PqO6/gJYarlrr2yuw3VNA/SLvw4Az1tfFtV8C/JVSrtaz06LLF2zrlFLKFfADih0sp7WeDcwG49baFahbCOEA1T1YwX7dAt8CD1u/6W8MNAfigC1Ac+vIAHeML72+1UbP9s/Ag9b1RwLLimxrpPX1g8BPurr3hAshqrzKDsW6Xyl1CugJfK+UWg2gtd4LLAb2AauAp7XWedaz0nHAamA/sNi6LMBLwAtKqXiMPtW51va5QJC1/QWgcPiWEEI4K1VdTwK7du2qt27danYZQohqRim1TWvdtbTl5AotIYSwAwlXIYSwAwlXIYSwAwlXIYSwAwlXIYSwAwlXIYSwAwlXIYSwAwlXIYSwAwlXIYSwAwlXIYSwAwlXIYSwAwlXIYSwAwlXIYSwAwlXIYSwAwlXIYSwAwlXIYSwAwlXIYSwAwlXIYSwAwlXIYSwAwlXIYSwAwlXIYSwAwlXIYSwg2p7a22l1EXguAN3WQu45MD9OVp1Pr7qfGwgx2drDbXWwaUtVG3D1dGUUlvLci/zqqo6H191PjaQ4zOLdAsIIYQdSLgKIYQdSLjazmyzC7Cz6nx81fnYQI7PFNLnKoQQdiBnrkIIYQcSrkIIYQcSriVQSj2klNqrlMpXSnW95rNXlFLxSqkDSql+Rdr7W9vilVIvF2lvrJTarJQ6pJRapJRyt7Z7WN/HWz9v5KjjK0opNUEpdVoptdP6uKvIZzY5VmdV0nFUBUqpY0qp3da/s63WtkCl1Brrn/8apVSAtV0ppWZYj3OXUqpzke2MtC5/SCk10qRjmaeUuqCU2lOkzWbHopTqYv2zireuq+x+UFpreRTzAFoDLYFfgK5F2tsAvwEeQGPgMGCxPg4DTQB36zJtrOssBh62vv4X8KT19VPAv6yvHwYWmXSsE4DxxbTb7Fid8XGj46gKD+AYUOuatmnAy9bXLwP/sL6+C1gJKKAHsNnaHggcsT4HWF8HmHAsvYHOwB57HAsQB/S0rrMSuNPexyRnriXQWu/XWh8o5qMBwJda6yyt9VEgHuhmfcRrrY9orbOBL4EB1v8h+wBLrOvPBwYW2dZ86+slwG0O+R+17Gx5rM6o2OMwuabKKvozde3P2gJt2AT4K6XqAf2ANVrrRK31ZWAN0N/RRWutfwUSr2m2ybFYP6uptd6ojaRdgAN+LiVcyy8UOFnk/SlrW0ntQUCS1jr3mvartmX9PNm6vBnGWX/Fmlfw6xe2PVZnVNJxVBUa+EEptU0pNcbaVkdrfRbA+lzb2l7ev0tnYKtjCbW+vrbdrlztvQNnppT6EahbzEevaq2XlbRaMW2a4v+j0jdY/kbbsrkbHSvwT2CSdd+TgPeAx25QX0WO1RlVtXqvFaW1PqOUqg2sUUr9foNlSzrWqvhnUN5jMeUYb+pw1VrfXoHVTgH1i7wPA85YXxfXfgnj1xZX6xld0eULtnVKKeUK+HH9r0Y2UdZjVUrNAb67pr4ClTlWZ3Sj43N6Wusz1ucLSqmlGN0c55VS9bTWZ62/Dl+wLl7SsZ4C/nRN+y92Lr2sbHUsp6yvr13erqRboPy+BR62ftPfGGiO0Vm+BWhu/bbcHeMLqm+tfTw/Aw9a1x8JLCuyrYJvNB8EfrIu71DWH9wC9wMF39ja8lidUbHHYXJNZaKU8lFK1Sh4DfTF+Hsr+jN17c/aCOs37T2AZOuv2quBvkqpAGt3UF9rmzOwybFYP0tRSvWwfi8wAkf8XDr6W8Gq8sAImVNAFnDe+pdU8NmrGN8yH6DIt44Y32IetH72apH2JhihFA/8D/Cwtnta38dbP29i0rF+BuwGdmH84Naz9bE666Ok43D2h/XP+TfrY29B7Rj93jHAIetzoLVdATOtx7mbq0fAPGb9+4oHHjXpeL4AzgI51n93o215LEBXjP98DgMfY7061Z4PufxVCCHsQLoFhBDCDiRchRDCDiRchRDCDiRchRDCDiRchRDCDiRchRDCDiRchRDCDv4/lGVvGZDOz4kAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Generate towers with x and y coordinates.\n",
    "# for tower i: x, y = towers[i][0], towers[i][1]\n",
    "towers = (np.random.rand(num_towers, 2)-0.5) * rx_square_side\n",
    "print('towers:\\n', towers)\n",
    "\n",
    "# location of transmitting device with tx[0] being x and tx[1] being y.\n",
    "tx = (np.random.rand(2)-0.5) * tx_square_side\n",
    "print('tx:', tx)\n",
    "\n",
    "# Distances from each tower to the transmitting device,\n",
    "# simply triangle hypotenuse.\n",
    "# distances[i] is distance from tower i to transmitter.\n",
    "distances = np.array([ ( (x[0]-tx[0])**2 + (x[1]-tx[1])**2 )**0.5\n",
    "                       for x in towers])\n",
    "print('distances:', distances)\n",
    "\n",
    "# Time at which each tower receives the transmission.\n",
    "rec_times = distances/v + t_0\n",
    "# Add noise to receive times\n",
    "rec_times += np.random.normal(loc=0, scale=rec_time_noise_stdd,\n",
    "                              size=num_towers)\n",
    "print('rec_times:', rec_times)\n",
    "\n",
    "# Get the loci.\n",
    "loci = get_loci(rec_times, towers, v, delta_d, max_d)\n",
    "\n",
    "# Plot towers and transmission location.\n",
    "fig, ax = plt.subplots(figsize=(5,5))\n",
    "max_width = max(tx_square_side, rx_square_side)/2\n",
    "ax.set_ylim((max_width*-1, max_width))\n",
    "ax.set_xlim((max_width*-1, max_width))\n",
    "for i in range(towers.shape[0]):\n",
    "    x = towers[i][0]\n",
    "    y = towers[i][1]\n",
    "    ax.scatter(x, y)\n",
    "    ax.annotate('Tower '+str(i), (x, y))\n",
    "ax.scatter(tx[0], tx[1])\n",
    "ax.annotate('Tx', (tx[0], tx[1]))\n",
    "\n",
    "# Iterate over every unique combination of towers and plot nifty stuff.\n",
    "for i in range(num_towers):\n",
    "    if(plot_trilateration_circles):\n",
    "        # Circle from tower i to tx site\n",
    "        circle1 = (towers[i][0], towers[i][1], distances[i])\n",
    "        circle = plt.Circle((circle1[0], circle1[1]),\n",
    "                            radius=circle1[2], fill=False)\n",
    "        ax.add_artist(circle)\n",
    "    for j in range(i+1, num_towers):\n",
    "        if(plot_lines_between_towers):\n",
    "            # Line between towers\n",
    "            ax.plot((towers[i][0], towers[j][0]),\n",
    "                    (towers[i][1], towers[j][1]))\n",
    "\n",
    "for locus in loci:\n",
    "    ax.plot(locus[0], locus[1])\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVcAAAEyCAYAAABOCJPxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xd0FFUbx/HvTU/oJdRA6L0TmoAgKFUFRRFRitKkg/Teq/QivYqKFUSkSJMmLfQqhB56b+nJff/YgTdigITs7uyG53POnt29O3PnGRJ+mb3TlNYaIYQQ1uVidgFCCJEUSbgKIYQNSLgKIYQNSLgKIYQNSLgKIYQNSLgKIYQNSLgKIYQNSLgKIYQNSLgKIYQNuJldgK2kT59e58iRw+wyBHDj4i3u3bxPnpI5zS4l3q6eu879Ww/xy5cZnxTeZpcjHMjevXtvaq19XzRdkg3XHDlyEBgYaHYZApjRdSFr5m/kt8DFZpcSL79MXMnMbotoOrghTQZ9aHY5wsEopc7HZzoZFhA25+LqQkx0jNllxEvgnweZ3WMxlRuU45MBDcwuRzgxCVdhcy6uLkQ7QbheCrrCiEYT8S+cjR4L2uPiIv89xMuT3x5hc86w5frofggD643BxdWFIct74p1cxllF4iTZMVfhOFxcFDrGccM1JiaG0U2mEHzyCmP+HEDmnBnNLkkkAbLlKmzOxdWFmBiNo147eOGApez8fS/tJn1GiTeKmF2OSCIkXIXNubhafs1iHHDrdcfvgXw/ahl1Wlbn3XY1zS5HJCESrsLmnoSrg427hj4KY1rHeeQsmp0O01qglDK7JJGEyJirsDlXBw3X70f+yvULN5mweSjuHu5mlyOSGNlyFTb3/2EBxxlzDT55mZ/GreCtplUoWrmg2eWIJEjCVdjck3CNija5EgutNVM7zMXD24NWYz41uxyRREm4Cpvz9PEEICwkwuRKLDYt3c6+9Yf5bPjHpMmY2uxyRBIl4Sps7vGFT0Luh5hcCTy8+4iZXy4kX0Bu3mlbw+xyRBIm4SpszielEa4PwkyuBOb3/Y57N+7TZWZrXF1dzS5HJGESrsLmHm+5hj4INbWO47tOsXLWOup3rEPeUrlMrUUkfRKuwua8U3gBEHLfvHCNjopm0hezSJclDc2GfmRaHeLVIeEqbM4npQ8AISZuuS6bsoozB8/TbvLncvFrYRcSrsLmfIwtV7OGBa5fvMmiQT9Qrm4pKr1X1pQaxKtHwlXYnLexpfjIpGGBrzvPR8doOkyVU1yF/Ui4Cpvz8vHExUWZsuX694o9bF++hyaDGpIpRwa7L1+8uiRchc0ppfBO4W33MdfQR2FM7zSfHEWy0aBrXbsuWwi5cIuwC5+U3jyy80kE3w77mesXbjJx6zDc3OVXXdiXbLkKu0ibKTV3rt612/KCT13hl4krqdG8KkUqFrDbcoV4zCrhqpSar5S6rpQ6EqstrVJqnVLqlPGcxmhXSqkpSqkgpdQhpVSpWPM0M6Y/pZRqFqu9tFLqsDHPFCV7JZxOer903Ay+bbflzeq+CA8vD1qMbGy3ZQoRm7W2XBcCtZ5q6w1s0FrnBTYY7wFqA3mNR2tgBljCGBgElAPKAoMeB7IxTetY8z29LOHgfLOm40bwLbssa8/aA+z8fS+f9G9A2kxpXjyDEDZglXDVWm8Bnt4sqQcsMl4vAurHal+sLXYCqZVSmYGawDqt9W2t9R1gHVDL+Cyl1nqHttyEaXGsvoSTSO+Xjkf3QuyyU2vx4B/JnCsj9TvVsfmyhHgWW465ZtRaXwEwnh8fB5MVuBhrumCj7XntwXG0/4dSqrVSKlApFXjjxg2rrISwDl+/tAA233oNOnCWE7tOUb9jbTw85e4Cwjxm7NCKa7xUv0T7fxu1nq21DtBaB/j6+iaiRGFt6f3SAXDTxuH6x6x1eHi581bTKjZdjhAvYstwvWZ8pcd4vm60BwPZYk3nB1x+QbtfHO3Cifga4XrDhju1Qh6EsuHbrVT56DVSpElus+UIER+2DNcVwOM9/s2A32K1NzWOGigP3DOGDdYCNZRSaYwdWTWAtcZnD5RS5Y2jBJrG6ks4iXRZLcMCttxy3fT9NkIfhvF267dstgwh4ssqR1Yrpb4HqgLplVLBWPb6jwZ+VEq1AC4AHxqTrwLqAEFACPAZgNb6tlJqGLDHmG6o1vrxZk5bLEckeAOrjYdwIh6e7qTOkMpm4aq1ZuWsdeQq5k/B8vlssgwhEsIq4aq1/vgZH1WPY1oNtH9GP/OB+XG0BwJFElOjMJ+vX1puXLJNuJ4MPE3Q/rN0nNZSLs4iHIKcoSXsJr1fOq6fv2mTvtfM34iXjyfVP6lkk/6FSCgJV2E3/gX9CD55mciISKv3fWDTEUpUL0KyVMms3rcQL0PCVdhN7hI5iIqM5vzR4BdPnAC3rtwh+OQVilUuZNV+hUgMCVdhN3lK5gQgaP9Zq/Z7ZOtxAIq+XtCq/QqRGBKuwm6y5MmEd3IvTh84Z9V+D205hlcyzyfhLYQjkHAVduPi4kKu4v4EHbDuluvhrccp9Fp+uWarcCgSrsKuchfPwZmD54mJibFKf/dvP+Ds4Qsy3iocjoSrsKs8JXMS8iCUK2euWaW/I9tOAFCsioSrcCwSrsKu/r9T65xV+ju85Tjunu7kL5PbKv0JYS0SrsKu/Atnw9XNldNWGnc9vPUYBcrlwcPLwyr9CWEtEq7Crjw83fEv5EeQFY4YCHkQyql9Z2W8VTgkCVdhd7lL5uDU3jNYLjPx8k4GniYmOobCleQGhMLxSLgKuytZrSh3r9/j1L4zieonIsxyGm2yVD7WKEsIq5JwFXZXtnZJlFLsWrkvUf24ull+fWOioq1RlhBWJeEq7C5V+pQUrJCPHSsDE9WPq5srANFR1jlmVghrknAVpihftzSn9p7h5uWXv+3L/8NVtlyF45FwFaYo/3YpAHav2v/Sfbi6S7gKxyXhKkyRo0h2MmRPz+5Ve1+6DxkWEI5MwlWYQilFuTql2LvuEBHhL3fx7Mc7tGTLVTgiCVdhmnJ1SxH2KJxDm4+91Pwy5iocmYSrME3xN4rg4eXOrj9ebmhAhgWEI5NwFabx8vGkZPWi7Ppj30udrSXDAsKRSbgKU5WrU4orZ65x9vCFBM/r7mG5OHZ4SLi1yxIi0SRchamqNHwNDy93Vkxfk+B502VNi08Kb6vfNkYIa5BwFaZKmS4F1T95nfVLtnD/9oMEzevq6kqB8nk5uuMfG1UnxMuTcBWmq9+xNuGhEayZtzHB8xaukJ9zhy8Q8iDUBpUJ8fIkXIXpchXzp1iVQqz4ei3R0QnbOVWwQj5iYjQndgfZqDohXo6Eq3AI9TvW4dr5G+z8PWGHZRUslxelFMf+lqEB4VgkXIVDeO3dADJkT8/yqasSNF/y1MnwL+zHMRl3FQ5GwlU4BFc3V95tV4sDm45y9vD5BM1buEJ+ju04abXbdQthDRKuwmHUblkNT28Plk9dnaD5ClbIx6N7IVw4fslGlQmRcBKuwmGkTJuCao0rs+Hbrdy/Ff/Dsgq/lh+AYztO2qo0IRJMwlU4lMeHZa1OwGFZWfNmJmW6FLJTSzgUCVfhUHIV86d41cIsn7qKiLCIeM2jlKLQa/k4tOVYou8oK4S1SLgKh/NJ/wbcvHQ7QWOv5euW5sqZa5w5lLCdYULYioSrcDglqxWlTO2SfD9qWbxPia34XllcXF3Y/OPfNq5OiPiRcBUOqeWoT3h0L4Slo5bFa/rUvqkoUa0Im3/aIUMDwiFIuAqHlKuYP282fZ3lU1dz7fyNeM1T5cPXuBx0Va6SJRyChKtwWM2HNgKlWDhwabymr1i/jGVo4KcdNq5MiBeTcBUOK0O29LzXqQ4blmzl9MFzL5w+VfqUlKxelC0//S1DA8J0Eq7CoTXqXZ/kqX2Y23tJvKav8mEFLp++RtD+szauTIjnk3AVDi1FmuQ07teAwLUH2bvu4Aunr1i/LK5urnLUgDBdkg3Xezfvc+/mfbPLEFbwbvtaZMqZgakd5hEe+vz7ZaVMl4JSbxVjw7dbiYyItFOFQvxXkg3Xa+dv0jBzK3rXGs6a+RsTfAsR4Tg8PN3pOqsNl05dYfHgn144/Xud6nDz0m02LNlqh+qEiJtKqgP/RQsV1V3q92Hzj39z5cw1XN1cKf5GYSq9V46K9cuQNlMas0sUCTSh1UzWLtjIxK3DKFQh/zOn01rTvkwvQh6EMe/YRFxdXe1YpUjqlFJ7tdYBL5zO1uGqlDoHPACigSitdYBSKi3wA5ADOAc01FrfUUopYDJQBwgBmmut9xn9NAP6G90O11ovet5yAwICdGBgIFprTu49w7ZfdrL1111cOnUFpRSFK+a3BO17ZcmUI4P1V1xY3aP7IXxRojvKxYWZ+7/CJ4X3M6fd+usuhn4wjr7fdeGNRhXtWKVI6hwtXAO01jdjtY0FbmutRyulegNptNa9lFJ1gI5YwrUcMFlrXc4I40AgANDAXqC01vrOs5b7OFxj01pz7uhFtv26i23LdnHmoOU89LylclLp/fJUer8c2QtkteLaC2s7sv0E3aoM5M2mVegxv/0zp4uJiaF1sW4oF8WsA+NwcUmyI2DCzuIbrmb9xtUDHm95LgLqx2pfrC12AqmVUpmBmsA6rfVtI1DXAbUSulClFDmLZKfJwA+ZtX8cC09OodWYT3HzcGNB/+9pUagLLQp3YUH/7zm174wcK+mAilQswMd93+fPhX8992QBFxcXPu7zPueOXGTHisBnTieErdhjy/UscAfLFucsrfVspdRdrXXqWNPc0VqnUUqtBEZrrbcZ7RuAXkBVwEtrPdxoHwCEaq3HPWu5cW25Ps+N4FtsX76b7ct2cWjzMWJiNJly+FLxvXJUer8chSrkk60fBxEVGUXXygMIPnmFWQfHkSFb+jini46K5vOCnUmeJjnTdo3CMuokROI40pZrRa11KaA20F4p9fpzpo3rt18/p/3fMyvVWikVqJQKvHEjfuejP+brl476HWrz1YbB/HBlDt3mtsW/cDZWTF9D18oD+NivDZPbzmbvuoNERUYlqG9hXW7ubvRe0omoyCi+aj7tmffOcnVzpVHv9zgZeJrAP198jKwQ1mTXowWUUoOBh0AroKrW+orxtf8vrXV+pdQs4/X3xvT/YNlqrWpM38Zo/9d0cUnoluuzPLofwu5V+9m2bBe7V+0j7FE4yVMno8K7AVR6rxylaxTD09sz0csRCbdm/kbGt5xBqzGf0rBHvTiniYyIpFnejmTInp6JW4bJ1qtINIfYoaWUSga4aK0fGK/XAUOB6sCtWDu00mqteyql6gId+P8OrSla67LGDq29QCmj631YdmjdftayrRWusYWHhrN33SG2/bqLHSsCeXj3EV7JPClTuySV6pelXN1SJEuVzKrLFM+mtWZYw/HsWBHIlB0jyVsqV5zTLZ+2mumd5jNu02CKVyls5ypFUuMo4ZoLeHxBTjfgO631CKVUOuBHIDtwAfhQa33bOBRrGpadVSHAZ1rrQKOvz4G+Rl8jtNYLnrdsW4RrbFGRURz86yjbft3F37/t4fbVu7i5u1KyelEqvVeOCvXKkCZDKpstX1jcv/WA1sW74ZPSh68Dx+Dl899vEeGh4TTJ1Z6cRbMz5s+BJlQpkhKHCFcz2TpcY4uJieH4zlNs+3UX25fv5sqZayilKFKpwJNjaTP6+9qlllfRvg2H6fXWUN75ogadvm4V5zQ/jVvB7J7fMGXHSAqWy2vnCkVSIuFqx3CNTWvNmUPn2b5sN9uW7eLs4QuA5VjaivXLUalBOfwL+tm9rqRuVvfF/Dzhd4b+1osK7/z39z70YSif5GhH4Yr5GfZbbxMqFEmFhKtJ4fq0S0FXngTt8Z2nAMutoMvXLUW5t0tTpFIB3D3cTa7S+UWER9KxfB9uXbrN7EPj4zy9ecmwn1k06Adm7BtLnhI5TahSJAUSrg4SrrHdvHSL7cv3sOuPvRzYeITIiCh8UnoTULME5d8uTdnaJUmVPqXZZTqt88cu0i6gF8WrFmb4yj7/OS754d1HfJKjLaXfKsbAn7qbVKVwdhKuDhiusYU+DGX/hiPs+D2Q3av2cfvqXZRSFKyQjwpvl6bc26XJUTibHDqUQCu+XsvUDnNpO7E573eu+5/PFw5cyrfDf2H2ofHkLJLdhAqFs5NwdfBwjS0mJoZT+86ya+Vedq4M5NQ+y1X0M/r7UrZ2ScrVLUXxN4rEuSdc/JvWmoH1x7Bn9QGGr+xDQI3i//r8/u0HfJqjHeXqlqLf911NqlI4MwlXJwrXp928dItdf+xj16p97F9/mLCQcDy83ClRrQjl6pSmbJ2SciWv53h0P4QvqwzkyulrTNgy9D/jq/P6fscPY5Yz58gE2bkoEkzC1YnDNbaIsAgObj7G7lX72L1qH5dPXwPAv5Af5eqUomzdUhR+LT9u7m4mV+pYbl66RacK/YiOjmHqjhFkyP7/Q+Hu3bzPpznbUbF+WXp/08nEKoUzknBNIuEam9aa4JOX2b1qP7tW7ePwlmNERUaTLJUPpWsUp0ytkpSpVYJ0meVC4ABnj1yga+UB+PqlY+LWYSRP/f+z5+b0/IafJ/zOvGOT8MuXxcQqhbORcE2C4fq0R/dD2Lf+sGWrdvV+bl+xXN42V3F/ytQsQZnaJV/5rdr9Gw/Tt/YIClcswMjV/fDwtBz2dufaXZrkas/rDSvQc0EHk6sUzkTC9RUI19gen7ywZ80BAtce4Mi2E0RHReOTwpuS1Ys82aqN/fX4VbF+yRbGNJ1KtcaV6LW445NDtGZ0XcjyaatZcGIyWXJnMrlK4SwkXF+xcH3ao/shHNh4hD2r97Nn7QGuX7DcCCJ7wayUqVmCgFolKfZ6QTy8PEyu1D6+G/krC/p/T6Ne9Wkx6hMAbl6+TdPcHaj+SWW6zW1rcoXCWUi4vuLhGpvWmgsnLj0J2sObjxEZEYWHlzvFqhSiTM2SlK5RjOwF/ZLscbVaayZ/MZs/5qyn09eteOeLGgBM6ziPlbPWMevAV/gXymZylcIZSLhKuD5T6KMwDm0+RuDaA+z98yAX/7kMQJqMqShetTDFqxaheNVC+OXLkqTCNjoqmkHvjWXP6v0M/6MvZWqW4M71e7Qs3JVsBbLI9V5FvEi4SrjG27XzN9i77hCHNh/lwKYj3Lps2TGWNnMailctRLHXC1P09YJkL5DV6cMn9GEoXSoN4Nr5G0zdOZJs+bOyet4GJrSaSe9vOlH9k8pmlygcnISrhOtL0VpzKegqh/46yoG/jnDwr2NPjkJI7ZuSwpUKUKxyIYpULkDu4jlwdXM1ueKEu3ruOh3K9iZF2uRM2TGSZKl86Fi+Lzcv3WbBicnPvWW3EBKuEq5W8Thsj2w9zqGtxzi85ThXz14HwCeFN4Ur5qdIpYIUe70g+crkeXKok6M7vPU4Pd8cQvE3ijBiZR/+CTxN59f6/WuHlxBxkXCVcLWZG8G3OLz1OEe2Hufw1uOcO3oRAHdPdwqUzUPB8vkoWD4vBcvnc+gTGlbN3cDE1jN5v3Nd2k5sztjPprHpu23MOTIRv7yZzS5POCgJVwlXu7l/6wFHtp3g0JZjHN1+gqD9Z4mKjAYgQ/b0lqAtZwncPCVzOtThX9M7z2f51NU0HdyQOq3e5PMCnSn6ekGG/97H7NKEg5JwlXA1TURYBEEHznFi5ymO7zrJ8Z2nuHbecqtzN3dX8pTMSYFyli3bfKVzkSVPpv9ce9VeoqOjGd9iBusWb6bpoIZ4JfNkds9vGP57b8rVLW1KTcKxSbhKuDqUW1fucGLXKY7vPMnxXac4uec0YSHhAHgn9yJXcX/ylMhJnlK5yFMyB/6F/Ox2h4bYAdu4XwO2/ryD6OgY5hye4DRjyMJ+JFwlXB1adFQ0Z49cIGj/OYL2nSHowFlOHzhH2CNL4Lp7uOFfOBt5SuQwAjcnOYpkI1lKH9vUEx3N+JYzWLdoM9UaV2bjd1tpMeoTGvWqb5PlCecl4Srh6nRiYmK4dOqKJXD3nyXowFmC9p3l/q0HT6ZJnzUt2Qv5kb1AVvwLZSN7waz4F/Kzyu1xoqOjmdBqJn8u/Au/fFm4eekWC/6ZQvosaRPdt0g6JFwlXJMErTU3gm8RtP8s548Gc/74RS4cv8TF45eeDCsApEqfwghdP/wL+eGXPwuZc2Uko3/6BA0vREdHM7HVLNYu3ISLqwtVP3qNPks622LVhJOScJVwTdJiYmK4cfEWF44Hc/5YMBeOX+L88WAuHg/mwZ1HT6ZTSpHeLy2Zc2UkU84MZM5pPBvv02ZK/Z+zzmJiYpjQciZrF24CYPLfIyhUPp9d1084rviG66t7oU/h1FxcXMjo70tGf1/K1Cr5pF1rzZ1rd7l06iqXT1/l6tnrXD13nStnrrH3z4NPTu19zNPbgwz+vqTPmpZ0WdKQLnMa0mVJS5laJbh67hoH/zrGmCZTWHhyqtOf+iteTGtNeGgEoQ/DCHsYRsiDUEIfhBLyIMx4Do13X7LlKl4p4aHhXDt/k6tnLYF79ex1rp2/zq3Ld4zH7SfH6Mbm4eVOljyZSJs5DanSpyBl2hSkTPf4kZwU6VKQMm3yJ+99UvpIGNuI1prI8EjCQsIJD4kgPDSCsEdhhD0KNx7/fx36MCzW+zDCQmK1PwwjNNbj8fsXZeJ6/bNsuQrxNE9vT7IXyEr2Alnj/FxrzYPbD7l1+TY3L9/h0qmrzPxyARFhkSRP5UPIvRCunL7Gg9sPeXj3UZx9ALi6uZIibXKSpfLBJ4UXXsm98EnhjXdyL7yS/f+195NnL7x8PHH3dMfd0w0PLw/cPd2M9+5PXnvEeu3i6mKX44O11kRHRRMTHUN0VDTRUY+fLa8jwyOJjIgiynjEfh8ZHhXrdSQRYZFEhEU8eY4Mi9UW+/NQS2hGhEYQFhJBeEi48bC0J3Sj0N3DDa9knngl88IrmeeTf/f0WdNafjbJvSw/m+ReeCc3fibJH//cLD8nnxTe+KTwIqP/z/FapoSrELEopZ5skeYs6k+ZmpAslTdjm03j1L6zTN05kpxF/QHL4WQP7jzk/q2H3L/1gAe3Lc9P3t96wKMHoU++Xt6+epfQB6GWLaUH8dtKig9LyCrL8+OHi8u/3qP1k2VpzZP3WlvCk1ivdYyOFaTRxMTY7tutp7cHHl7uuHtZnj283PH09sDd0x1PH09SpEmOp48Hnj6eeHoZzz4eeHpbnr18PPHwtrRb/nB5/jtEjddmXGBIhgWEeAGtNR3L9+Xk3tOkyZCKiVuHWeW2MFpry9dUI3zDHoVbtvrCja2/8EgiYr2OfOp1THQM0dHRT8LQ8j7m/+FotMVEx6CU5Q8HxlBF7PePX1veW8azXd1ccHVzffJw+dd7y2sXV8uzu6cb7h5uuHkYW9rGazcPS7u7p/uT148D1N3LA3cPN6ccOpGjBSRchRWd2H2KjuX74uHlTvI0yRm9pt+TLVjxaolvuJpzQrcQTqZA2bxUaVgBlELHxPBllUEc2/GP2WUJBybhKkQ8fT6iMTpGk6t4DlKkS06vt4axZ+0Bs8sSDkrCVYh4ypI7E58N/5i9fx7k/U51yJI3EwPfHc1fP2w3uzThgCRchUiA97vUoUilAiwc+AN9vulEgXJ5Gdl4MitnrTO7NOFgJFyFSABXV1d6LGhPdFQ0c3ovYeTqfpStU5LJbWezfNpqs8sTDkTCVYgEypI7E82GfMTuVfsJXHuQwb/24LV6ZZjeaT6r5qw3uzzhICRchXgJ73WqQ67i/nzdeT4RYZH0W9qVMrVLMumL2axbvNns8oQDkHAV4iW4urnSeUZrbl2+w+LBP+Lh6c6gn7tRoloRxn0+XXZyCQlXIV5WofL5qNOyOsumrOL0wXN4ensyZHlPClcswKhPp7Bt2S6zSxQmknAVIhE+H9WYlGmTM7ntbGJiYvBO5sXwlX3IXyY3IxpN5O8Ve8wuUZhEwlWIREiZNgWtv2rK8Z2nWD13AwA+KbwZuaofuUvkYEiDcayet8HkKoUZJFyFSKQ3m7xOsSqFmNfnW+5cvwdA8tTJGLt+EKXeLMqEVjP5ZuhPVrkClnAeEq5CJJJSik5ftyL0YRhzen7zpN0nhTfDVvTmrWZVWDz4Rya1mUV01H8vxC2SJglXIazAv6AfH3Z/l3WLN3Pwr6NP2t3c3egxvz2N+77PqrkbGPz+V4Q+CjOxUmEvEq5CWEnjfg3IlDMDk9vNITIi8km7UorPhn9Mp+kt2b1qHz2rD+HujXsmVirsQcJVCCvx8vGkw9QWXDxxiV8nrfrP5++0rcnAn7tz5tB5ulUdxO2rd+LoRSQVThOuSqlaSql/lFJBSqneZtcjRFzK1SlF+bdL892IX57s3IqtYv2yjFzdj+sXbtK92hBuXZGATaqcIlyVUq7AdKA2UAj4WClVyNyqhIhb66+aEB4awaKBP8T5efEqhRm5qh83Lt6k+xuDuHnplp0rFPbgFOEKlAWCtNZntNYRwFKgnsk1CRGnbPmz8m67mqyeu56zh8/HOU3RygUZtaY/ty7fodsbg7kRLAGb1DhLuGYFLsZ6H2y0/YtSqrVSKlApFXjjxg27FSfE0z4d+AHJUidjcrs5REfHffhVkYoFGLWmP3ev3aP7G4O4fvGmnasUtuQs4RrXLSL/c0S21nq21jpAax3g6+trh7KEiFvKtCloO6E5R7f/w/Ipz77Oa+HX8jNqbX/u3rhPtyoDuXDikh2rFLbkLOEaDGSL9d4PuGxSLULEy5tNXqfCuwHM7/cdF/95dmgWKp+PsesHERYSQZeK/Ti05ZgdqxS24izhugfIq5TKqZTyABoBK0yuSYjnUkrRZWZrPH08Gdt8+jOHBwDyB+Rmyt8jSJ0hFb1rDGPTUrlkobNzinDVWkcBHYC1wHHgR6310efPJYT50mZKQ8dpLTmx6xQ/jfv9udNmzpWRSduHG/flmsTSMcvlegROzCnCFUBrvUprnU9rnVtrPcLseoSIr6ofvUblBuVYPOgHzh29+NxpU6ZNwei1/anaqCLz+nzL5LZz5HoETsppwlUIZ/X4wi7JUvkwtvk0oiKjnju9h5cHfZZ0olHYTbvdAAAaf0lEQVSv+vwxex0D648h9GGonaoV1iLhKoQdpPZNRaevW3Fq7xl+nrDyhdO7uLjQYtQndJ7RmsA1B+hVYxgP7z6yQ6XCWiRchbCTyg3KU/G9siwZ+hNXzl6L1zxvt3mLgT9359TeM/SoPoR7N+/buEphLRKuQthR+8mf4+LqwrSO8+K9s6pi/bIMWd6LC8eD6V5tMHeu3bVxlcIaJFyFsCNfv3Q0G/IRu1ftZ+svO+M9X9naJRm+sg9Xz1ynW1W5HoEzkHAVws7qd6xN7hI5+LrLAh7dD4n3fCWrFWXUWsv1CL6sMojgU1dsWKVILAlXIezM1c2VLjNbc/vKXRYOWJqgeYtULMDY9QMJuR9C59f6cWTbcRtVKRJLwlUIExQom5e3v6jBiulrOLn3dILmzV8mD5P/HkHKdMnp+eZQOZvLQUm4CmGSFiM/JnWGVEz6YvZzT42NS9Y8mZm8fcSTs7m+G/mrnM3lYCRchTBJslTJaDuxOaf2nmHF9LUJnj9luhSM/nMA1RpXYkH/75nQauYLT1AQ9iPhKoSJqjR8jdI1irNwwNKXOgLAw9Od3t904pP+DVgzfyN964yUkw0chISrECZSStFpekuiIqOY3nnBS321V0rRfGgjus9vx6HNx+hSqT/XzsvF4s0m4SqEybLkzkSTgR+y7dddido5VbP5G4xa04+bl27TuWI/zh8PtmKVIqEkXIVwAB/2eJcCZfMwvdP8RJ2BVbJaUSZuGUpMdAzdqgwkaP9ZK1YpEkLCVQgH4OrqSvf57Qh9EMrUjvMS1VfOov5M2DIMD28PulcbzLEd/1ipSpEQEq5COAj/QtloMqghW3/eyZafdySqL7+8mZm0dRipfFPSq8YwAv88aKUqRXxJuArhQBr2eJe8pXMxtf3cRF8BK0N2XyZuGUqWPJnoV3ckq+ast1KVIj4kXIVwIK5urvSY346Hdx8xvfP8RPeXNlMaJm4ZRukaxZnYZhZzei0hJibGCpWKF5FwFcLB5CzqT+N+Ddj0/Xa2L9+d6P58Ungz7LdevPNFDX786jeGfzSB8NBwK1QqnkfCVQgH9HGf98hdIgeT287m/u0Hie7P1c2VjtNb8sX4Zmz7dTfdqw2R68LamISrEA7Izd2N7vPbcf/WQ2Z0XWiVPpVSNOj6NoN+6c7ZQ+fpVKEv5489/4aJ4uVJuArhoPKUyEmj3vVZ/80Wdv2x12r9VqxflvGbhxIRFknniv3Zt+Gw1foW/yfhKoQD+6R/A3IUycbENrOses2A/AG5mbJjJL7Z0tG39gjWzN9otb6FhYSrEA7M3cOd7vPbc+faPeb2WmLVvjP6+zJp6zBKVCvC+JYz+H7UMqv2/6qTcBXCweUPyM17HWvzx5z1HNt50qp9J0uVjOG/96Za40rM7/cdc3svkevCWomEqxBOoOmQj0ifNS2T284mOiphF9Z+ETd3N3ot7sjbbd7ih7G/MbH1LLkurBVIuArhBHxSeNNu0mecOXie5VNXW71/FxcXOn3disZ932f1vA30f2c0j+7JdWETQ8JVCCdR6f1ylKldkkWDfuBGsPVvra2U4rPhH/PlnC84sPEIXSoNkOvCJoKEqxBOQilFx6ktiI6KZkbXBTZbTu0W1Rm5uh83gm/RsXwf/tkTZLNlJWUSrkI4kcy5MtK4XwO2/rKLrb/ustlySlUvyuTtw/H09qBb1UFsW2a7ZSVVEq5COJmGPd4lb6mcTGozi1tX7thsOf6FsjFlx0hyFvNn6Afj+Wn873IkQQJIuArhZNw93On1TSfCHoUxvuUMmwZemoypGbdxEJUalGN2j8VMaTfH6kcrJFUSrkI4If+CfrQa24Q9q/ezcuafNl2Wp7cn/Zd25aOe9Vg5ax0D3h3No/shNl1mUiDhKoSTqte+FgE1izOr+2Iu/nPJpstycXGh5ehP6Tr7C/auO0TXygO4fkGOJHgeCVchnJRSim7z2uHh7cFXn00nOtr2X9frtKzOqNX9uHb+Bl0qDeDCCduGujOTcBXCiaXPkpb2kz/n+M5T/DZtjV2WWerNYkzcMoyoyCi6Vh7Ake0n7LJcZyPhKoSTq9a4EmXrlGRBv++5cvaaXZaZq5g/E7cOI0Xa5PR8cyiblm63y3KdiYSrEE5OKUXnGa1xcXVhUptZdjtcKmuezEz5ewQFyuZhZONJfDviFzlUKxYJVyGSgAzZ0tNyzKfsW3+YtQs22W25KdOlYPSfA6j+SWUWDljKuBZfExkRabflOzIJVyGSiLqt36To6wWZ2W0RNy/ftttyPTzd6bW4I00GfsifC/+ib+0RPLjz0G7Ld1QSrkIkES4uLnw5py2R4ZFMbT/Xrl/RlVI0HdyQnos6cGTbCTpX7M+VM/YZ/3VUEq5CJCF+eTPTbMhH/P3bHrb8vNPuy3+rSRVG/zmAu9fu0qlCX47t+MfuNTgKCVchkpgGXd8mb+lcTOs4j/u3En9b7oQqXqUwU3aMxCeVD92rDWHzj3/bvQZHIOEqRBLj6uZK93nteHD7ITO+XGhKDX75sjDl7xHkC8jF8EYT+X7UslfuSAIJVyGSoFzF/J/clnv36v2m1JAqfUrGrhvIGx9XZH6/75jQcsYrdfsYm4WrUmqwUuqSUuqA8agT67M+SqkgpdQ/SqmasdprGW1BSqnesdpzKqV2KaVOKaV+UEp52KpuIZKKxv0a4F/Ij8ltZxMWEm5KDR5eHvRZ0plP+jdgzYJN9Kk9wqq3CHdktt5ynai1LmE8VgEopQoBjYDCQC3ga6WUq1LKFZgO1AYKAR8b0wKMMfrKC9wBWti4biGcnoenO51ntOb6hZssGviDaXUopWg+tBE9F3bgyNbjdK7Yz25nkpnJjGGBesBSrXW41vosEASUNR5BWuszWusIYClQTymlgGrAz8b8i4D6JtQthNMpWrkgdVu/xS8TV3Jk23FTa3mraRVGrx3A7St36VS+r9VvE+5obB2uHZRSh5RS85VSaYy2rMDFWNMEG23Pak8H3NVaRz3V/h9KqdZKqUClVOCNG3I5NCEAWn/VhIw5fBnbfDqhD0NNraV41cJM2TEC7xTe9Kg2mM0/7TC1HltKVLgqpdYrpY7E8agHzAByAyWAK8D4x7PF0ZV+ifb/Nmo9W2sdoLUO8PX1TfD6CJEU+aTwpvv8dlw9e505vb41uxyy5c/KlB0jyFMqF8M/msDS0UnzSIJEhavW+k2tdZE4Hr9pra9praO11jHAHCxf+8Gy5ZktVjd+wOXntN8EUiul3J5qF0LEU/EqhXm/cx1+n7GWvesOml0OqX1T8dX6gVT96DXm9f2OCa1mJrkjCWx5tEDmWG/fA44Yr1cAjZRSnkqpnEBeYDewB8hrHBnggWWn1wpt+ZO2CfjAmL8Z8Jut6hYiqfpsxMdky5+F8S1m8Oie+XvsPbw86PNtZxr3fZ818zfSt87IJHUkgS3HXMcqpQ4rpQ4BbwBdAbTWR4EfgWPAGqC9sYUbBXQA1gLHgR+NaQF6AV8qpYKwjMHOs2HdQiRJnt6e9FzUgVuXb/N114VmlwNYrofw2fCP6T6/HYc2H0tSRxKopDjWARAQEKADAwPNLkMIh7Og//d8N/JXhizvyWvvljG7nCcObDrCkAbjcPNwY+hvvShYLq/ZJcVJKbVXax3wounkDC0hXjGfDvyAXMX8mdRmlinXHniWEm8UYfLfI/BK5kn3Nwbx1w/OfXcDCVchXjHuHu70XNSBB7cfMrXDXLPL+ZfsBbIydedI8pTKxYiPJzG39xK73HjRFiRchXgF5S6eg08HfshfP/zN1l93mV3Ov6T2TcW4jYOo2/otfhj7G0M/GE/oozCzy0owCVchXlEf9axHruL+TOs4z+H20rt7uNNlZmvaT/6cnb8H0q3qIG5duWN2WQki4SrEK8rN3Y0v57Tl7rW7zO21xOxy4lS/Y20GL+vJxeOX6FShL2ePXDC7pHiTcBXiFZY/IDfvda7LH3PWc2jLMbPLiVOFdwKYsGUoUZHRdKnUn8A/zT8JIj4kXIV4xTUb+hGZcvgysfVMIsIizC4nTnlL5WLqzpFk9PelX92RrJqz3uySXkjCVYhXnHcyLzrPbEPwySt8N+JXs8t5pgzZ0jNx6zBKv1WMiW1mMbf3EmJiYswu65kkXIUQBNQozptNXmfpmOWcPXze7HKeKVlKH4at6M3bbSxHEgxvNJHwUHMuBP4iEq5CCAC+GN+M5Kl9GNfCsW/H4urmSqevW9H6q6Zs+2UXPaoP4c71e2aX9R8SrkIIwHLPq84zWnMy8DTfDPnJ7HKeSynFh93eYcBP3Thz8Dydyvfh9MFzZpf1LxKuQognKjcoT83mb7B09DLT71wQH5XfL8f4v4ZYjiSo2N+hbuMt4SqE+Jd2kz8jY44MjG4y1SEuTfgi+cvkYfqe0eQqkYPhjSYyr+93DrGjS8JVCPEvPim86f1NR25cvMnMbovNLide0mZKw1cbBlGnZXWWjl7GsIYTTLvj7WMSrkKI/yhUIT8Ne9RjzfyN7Fl7wOxy4sXD050us9rwxfhmbF+2mx7VBnPn2l3T6pFwFULEqcmgD8leMCsTW810iuEBsOzoatD1bQb+3I2zhy/QqUJfzh8PNqUWCVchRJw8vDzoPr89ty7fZlb3b8wuJ0EqvVeOcZsGExYSQZeK/Tmw6ciLZ7IyCVchxDMVLJeXD7q9y+p5G5zmnP7HCpTNy9SdI0mbOTV9ag1n3eLNdl2+hKsQ4rmaDWlItgJZmdBqBo/uh5hdToJkypGBydtHUKRyQcY2n8biwT/a7TbeEq5CiOeyDA+049al28zp4VzDAwDJUydj5Kq+1GhelW+G/sTY5tOICI+0+XIlXIUQL1SofD4++PId/piznn3rD5ldToK5e7jTfV47mg9txPpvttD9DdtffFvCVQgRL02HNCRb/iyMb+l8wwNgOZLgk/4NGPhTN84eukD7Mr34Z0+QzZYn4SqEiBdPb0+6z2/HjYu3mNPTMe9cEB+VG5Rn8t8jcPdwo+vrA1m/ZItNliPhKoSIt0IV8tOg69v8MXudUw4PPJarmD/Tdo+mUIV8jGk6ldk9Flv9LrMSrkKIBGk+7CP88mVmYptZpp9imhip0qdk9Nr+vNuuJj+N/50hDcZZ9S6zEq5CiATx9Paky8w2XD17nSVDHfvShC/i5u5Gx2ktaT/lc3at3Ev3NwZz+6p1dnRJuAohEqx41cLU+uwNfhr/u013CtlL/Q6Wu8xeOBZMpwr9OH/sYqL7lHAVQryUNuObkTZTasY0m+awt1pJiArvBDB+8xAiwyPpXLE/+zceTlR/Eq5CiJeSPHUyus9vx8UTl5jf93uzy7GKfKVzM3XnSHz90tG39ohEnTIr4SqEeGml3ypOvfa1+HXyH6ZcHMUWMmT3ZdK2YRSrUoixzacxv9/LXXxbwlUIkSgtx3yKX77MjG0+zWkuTfgiyVIlY8QffanTsjrfj1rG4Pe/IuRBaIL6kHAVQiSKl48nPRd15Nal20zvssDscqzGzd2NLrPaWI4k+GMfnV/rx+XTV+M9v4SrECLRCpbLy8d93mfdos1sW7bL7HKsRilF/Q61GbWmP7cu36ZDuT7xnlfCVQhhFZ8MaECekjmZ/MVs7ly/Z3Y5VlWqelGm7R5N2kyp4z2PhKsQwircPdzptbgjj+6HMqnNLLtdN9VesuTOxOS/R8R7eglXIYTV5Cicjc+GNeLv3/aw8bttZpdjdclS+sR7WglXIYRVvd+1LgXL5+XrLguS3PBAQki4CiGsytXVlW5z2xL6IJTpneaZXY5pJFyFEFbnXygbjfs1YPOPO9i+fLfZ5ZhCwlUIYRONetcnVzF/prSfy8O7SePkgoSQcBVC2ISbuxvd5rXl7vV7zOq2yOxy7E7CVQhhM/lK5+bDbu+wZsEm9q47aHY5diXhKoSwqSaDPrTcuaD1LEIfJuz8fGcm4SqEsClPb0+6zW3L9Qs3mdH11RkekHAVQthckUoF+ahnPVbP25Ckrj3wPIkKV6XUh0qpo0qpGKVUwFOf9VFKBSml/lFK1YzVXstoC1JK9Y7VnlMptUspdUop9YNSysNo9zTeBxmf50hMzUIIczQd0pC8pXIyodVMbl6+bXY5NpfYLdcjwPvAv278rZQqBDQCCgO1gK+VUq5KKVdgOlAbKAR8bEwLMAaYqLXOC9wBWhjtLYA7Wus8wERjOiGEk3H3cKfPt52JCI1gYuuZSe7aA09LVLhqrY9rrf+J46N6wFKtdbjW+iwQBJQ1HkFa6zNa6whgKVBPKaWAasDPxvyLgPqx+no8UPMzUN2YXgjhZLLlz0rL0Z+ye9V+1i7YZHY5NmWrMdesQOzbJwYbbc9qTwfc1VpHPdX+r76Mz+8Z0wshnNC77WtSrEohZny5kOsXbphdjs28MFyVUuuVUkfieNR73mxxtOmXaH9eX3HV2lopFaiUCrxxI+n+0IRwZi4uLnSf146Y6BjGt0q6wwMvDFet9Zta6yJxPH57zmzBQLZY7/2Ay89pvwmkVkq5PdX+r76Mz1MBcY6Ga61na60DtNYBvr6+L1o1IYRJMufKSOuxTdi37hCr5qw3uxybsNWwwAqgkbGnPyeQF9gN7AHyGkcGeGDZ6bVCW/50bQI+MOZvBvwWq69mxusPgI06qf6pE+IVUrfNW5SsXpRZ3Rdz9dx1s8uxusQeivWeUioYqAD8oZRaC6C1Pgr8CBwD1gDttdbRxphpB2AtcBz40ZgWoBfwpVIqCMuY6uNrlc0D0hntXwJPDt8SQjgvFxcXus1ti1KK8S1nvNTtqx2ZSqobgQEBATowMNDsMoQQL7Bq7gYmtp5Jx2ktebddzRfPYDKl1F6tdcCLppMztIQQpqrdohoBNYszp+c3Cbp1taOTcBVCmEopxZdz2uLq7sq4z79OMsMDEq5CCNP5+qWj7cTmHN56nOVTVptdjlVIuAohHEKNZlUpV7cU8/t9x5Uz18wuJ9EkXIUQDkEpRecZrXF1c2VS29lOf3KBhKsQwmH4+qWjxahP2LfuEH8u+svschJFwlUI4VDe/uItilYuyNddFjj1tQckXIUQDsXFxYUeC9ujYzRffTbdaY8ekHAVQjiczDkz0nZicw5sOsqyyavMLuelSLgKIRxSrc+rUf6d0szr+x3njl588QwORsJVCOGQlFJ8OfsLkqX0ZkzTqURGRJpdUoJIuAohHFaajKnpMqsNQfvPsmTozy+ewYFIuAohHFrF+mWp0bwqS0cv49iOuO4q5ZgkXIUQDq/dpM/wzZaeMc2mEfoozOxy4kXCVQjh8JKl9KHHwvZcOX2N+X2+M7uceJFwFUI4heJVClOvfS2WT1vN4a3HzS7nhSRchRBO4/NRjcmUMwPjWnxNWEi42eU8l4SrEMJpeCfzotvctlwOusrCAUvNLue5JFyFEE6lxBtFeLvNW/w66Q+HPnpAwlUI4XRajW2Cb7Z0jGsxg4iwCLPLiZOEqxDC6fik8Kbr7C+4eOISiwf/aHY5cZJwFUI4pYAaxan1eTV+GreCf/YEmV3OfyTZW2srpW4A5+24yPTATTsuz96S8vol5XUDWT9r89da+75ooiQbrvamlAqMz73MnVVSXr+kvG4g62cWGRYQQggbkHAVQggbkHC1ntlmF2BjSXn9kvK6gayfKWTMVQghbEC2XIUQwgYkXIUQwgYkXJ9BKfWhUuqoUipGKRXw1Gd9lFJBSql/lFI1Y7XXMtqClFK9Y7XnVErtUkqdUkr9oJTyMNo9jfdBxuc57LV+sSmlBiulLimlDhiPOrE+s8q6OqpnrYczUEqdU0odNn5mgUZbWqXUOuPff51SKo3RrpRSU4z1PKSUKhWrn2bG9KeUUs1MWpf5SqnrSqkjsdqsti5KqdLGv1WQMa+y+UppreURxwMoCOQH/gICYrUXAg4CnkBO4DTgajxOA7kAD2OaQsY8PwKNjNczgbbG63bATON1I+AHk9Z1MNA9jnarrasjPp63Hs7wAM4B6Z9qGwv0Nl73BsYYr+sAqwEFlAd2Ge1pgTPGcxrjdRoT1uV1oBRwxBbrAuwGKhjzrAZq23qdZMv1GbTWx7XWcV1ypx6wVGsdrrU+CwQBZY1HkNb6jNY6AlgK1DP+QlYDHt9dbRFQP1Zfi4zXPwPV7fIXNf6sua6OKM71MLmmxIr9O/X079pibbETSK2UygzUBNZprW9rre8A64Ba9i5aa70FuP1Us1XWxfgspdZ6h7Yk7WLs8Hsp4ZpwWYHYN1EPNtqe1Z4OuKu1jnqq/V99GZ/fM6Y3QwfjK9b8x1+/sO66OqJnrYez0MCfSqm9SqnWRltGrfUVAOM5g9Ge0J+lI7DWumQ1Xj/dblNutl6AI1NKrQcyxfFRP631b8+aLY42Tdx/qPRzpn9eX1b3vHUFZgDDjGUPA8YDnz+nvpdZV0fkbPU+raLW+rJSKgOwTil14jnTPmtdnfHfIKHrYso6vtLhqrV+8yVmCwayxXrvB1w2XsfVfhPL1xY3Y4su9vSP+wpWSrkBqfjvVyOriO+6KqXmACufqu+xxKyrI3re+jk8rfVl4/m6UmoZlmGOa0qpzFrrK8bX4evG5M9a12Cg6lPtf9m49Piy1roEG6+fnt6mZFgg4VYAjYw9/TmBvFgGy/cAeY295R5YdlCtMMZ4NgEfGPM3A36L1dfjPZofABuN6e3K+MV97D3g8R5ba66rI4pzPUyuKV6UUsmUUikevwZqYPm5xf6devp3ramxp708cM/4qr0WqKGUSmMMB9Uw2hyBVdbF+OyBUqq8sV+gKfb4vbT3XkFneWAJmWAgHLhm/JAef9YPy17mf4i11xHLXsyTxmf9YrXnwhJKQcBPgKfR7mW8DzI+z2XSun4DHAYOYfnFzWztdXXUx7PWw9Efxr/zQeNx9HHtWMa9NwCnjOe0RrsCphvreZh/HwHzufHzCgI+M2l9vgeuAJHG/7sW1lwXIADLH5/TwDSMs1Nt+ZDTX4UQwgZkWEAIIWxAwlUIIWxAwlUIIWxAwlUIIWxAwlUIIWxAwlUIIWxAwlUIIWzgf9j4aiFQHSikAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot the same hyperbola using the derived expressions.\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(5,5))\n",
    "max_width = max(tx_square_side, rx_square_side)/2\n",
    "ax.set_ylim((max_width*-1, max_width))\n",
    "ax.set_xlim((max_width*-1, max_width))\n",
    "\n",
    "c = np.argmin(rec_times)  # Tower that received first.\n",
    "p_c = towers[c]\n",
    "t_c = rec_times[c]\n",
    "\n",
    "x = np.linspace(towers[i][0] - 50000, towers[i][0] + 50000, 100)\n",
    "y = np.linspace(towers[i][1] - 50000, towers[i][1] + 50000, 100)\n",
    "x, y = np.meshgrid(x, y)\n",
    "\n",
    "for i in range(num_towers):\n",
    "    if i == c:\n",
    "        continue\n",
    "        \n",
    "    p_i = towers[i]\n",
    "    t_i = rec_times[i]\n",
    "    \n",
    "    plt.contour(\n",
    "        x, y,\n",
    "        (\n",
    "           np.sqrt((x-p_c[0])**2 + (y-p_c[1])**2) \n",
    "         - np.sqrt((x-p_i[0])**2 + (y-p_i[1])**2) \n",
    "         + v*(t_i - t_c)\n",
    "        ),\n",
    "        [0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Actual emitter location:    (731.7, 2184.2) \n",
      "Calculated emitter locaion: (711.6, 2129.4)\n",
      "Error in metres:            58.4\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVcAAAEyCAYAAABOCJPxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xd8FHX+x/HXdzebCklIo4YmoYMgoQhnRSl6p1KVIigg3il6/tSzI0VQz1PxOD1FRQFRAbEhJ02KSCcg0kuAAIFAQkICIaR/f3/MBAMGkpDdnd3weT4e+9jd77TPUN6ZfOc7M0prjRBCCOeyWV2AEEJURhKuQgjhAhKuQgjhAhKuQgjhAhKuQgjhAhKuQgjhAhKuQgjhAhKuQgjhAhKuQgjhAj5WF+AqERERun79+laXIQAyjkB2BlRvaXUlZVKgCziQfgCN5prQa7Aru9UlCQ+yadOmk1rryNLmq7ThWr9+feLi4qwuQwD88HfYswCe9vy/j0JdyCNLH0ElKWb0nEHLCO/4gSDcRyl1qCzzSbeAEMVM2TqF1UdX81yH5yRYRYVIuAphWntsLe9veZ+/NPwL/Rr3s7oc4eUkXIUAUs+l8sKqF2gQ0oCXOr2EUsrqkoSXq7R9rkKUVaEu5KXVL3E65zQf3PYBgY5Aq0sSlYAcuQo3UKALrS7ikmbunMmqo6t4uv3TNAlrYnU5opKQcBWuZ3dAQZ7VVZRoR+oOJm2exK3Rt3Jfk/usLkdUIhKuwvXsvlCQa3UVf5CVl8WzK58l3D+c8V3GSz+rcCrpcxWu5+MH+TlWV/EH/4r7F4dPH2Zq96mE+IVYXY6oZOTIVbie3Rd0ARQWWF3Jeb8k/sLcvXMZ2mIo7Wu0t7ocUQlJuArXs/sa7x7SNZCRk8GYNWNoFNqIUW1HWV2OqKSkW0C4no+f8Z6fA44Aa2sBJq6fyKnsU7zb9V387H5WlyMqKTlyFa7nQUeuSw4tYcHBBTx87cM0D29udTmiEpNwFa7n42+8552ztIzTuad5df2rNAtrxohWIyytRVR+0i0gXM+vqvGem2lpGf/e9G/SstN4t+u7+Njkn75wLTlyFa5XFK7Zpy0r4dfkX5mzdw4Dmw6kRXgLy+oQVw8JV+F6/uYY0pwzlmw+ryCPcWvGUSOoBo+1fcySGsTVR343Eq5XdOSaY82R67Qd09ifsZ//3PofuSmLcBs5chWu5xdsvFsQrifOnuCjbR/RtW5Xbo6+2e3bF1cvCVfheuePXN3fLTD518nkF+bzVOxTbt+2uLpJuArX8w0CmwOy0ty62R2pO5i3fx6Dmw0mumq0W7cthISrcD2lICgCsk66bZNaa97Y8AZh/mE81Poht21XiCISrsI9giLhrPvCdWXiSjYnb+bRNo9S1beq27YrRBGnhKtS6hOlVLJSanuxtjCl1BKl1D7zvZrZrpRSk5VS8UqprUqp64otM9Scf59Samix9nZKqW3mMpOV3HjT+wRFwtkUt2xKa817W94jumo0vWJ6uWWbQlzMWUeu04AeF7U9ByzVWscAS83vAD2BGPM1EngfjDAGxgAdgQ7AmKJANucZWWy5i7clPJ0bw3Xp4aXsStvFX6/9Kw6bwy3bFOJiTglXrfVK4OKzFXcD083P04F7irXP0IZ1QKhSqibQHViitU7TWp8ClgA9zGnBWuu1WmsNzCi2LuEtgiIgMwW0dulmtNa8/9v71A+uzx0N7nDptoS4HFf2uVbXWicBmO9RZntt4Eix+RLNtsu1J5bQ/gdKqZFKqTilVFxKinuOkkQZBdeG/HMuHzGwNmkte0/tZVjLYXL/AGEpK05oldRfqq+g/Y+NWn+otY7VWsdGRkZWoEThdKHmUKiMwy7dzIydMwj3D+fOhne6dDtClMaV4XrC/JUe8z3ZbE8Eig86rAMcK6W9TgntwpuEmH+16UcuP18FxJ+KZ/XR1QxoOgDfonvICmERV4brPKDojP9Q4Pti7UPMUQOdgAyz22AR0E0pVc08kdUNWGROO6OU6mSOEhhSbF3CW4TWNd4zXBeun+/+HD+7H/2b9HfZNoQoK6d0SimlvgRuBiKUUokYZ/1fB+YopYYDh4F+5uw/AncA8UAW8CCA1jpNKfUKsNGcb7zWuqiD7m8YIxICgAXmS3iTgGrgCHLZkWtWXhY/HviR7vW7U82/WukLCOFiTglXrfWAS0zqWsK8Gnj0Euv5BPikhPY4oGVFahQWU8rod013TZ/r4kOLycrPok9MH5esX4jykiu0hPuENYS0/S5Z9YKDC6hTpQ5to9q6ZP1ClJeEq3CfiMaQuh8K8p262lPZp1iftJ4eDXogF+8JTyHhKtwnojEU5sGpBKeudvmR5RToArrV6+bU9QpRERKuwn0imxjvJ/c4dbVrjq0hKiCKpmFNnbpeISpCwlW4T0Rj4z1lt9NWWagLWZ+0nk61OkmXgPAoEq7CffyDIbQeJG112ip3p+0mPSedTjU7OW2dQjiDhKtwr1pt4divTlvd5hObAehYs6PT1imEM0i4Cveq1QbSDzntBi7JWcn42nyJDJB7SQjPIuEq3KuWOQ41aYtTVpeanUpYQJj0twqPI+Eq3Kvmtca7k7oG0rLTCPcPd8q6hHAmCVfhXgHVjCu1jm52yupSz6US5h/mlHUJ4UwSrsL9ojvBoTVQWFjhVaVlp0m4Co8k4Srcr34XOJdW4fGuWmsjXAMkXIXnkXAV7levi/F+aHWFVpOZl0leYZ70uQqPJOEq3K9afePJBAdWVGg1qedSAaRbQHgkCVfhfkpBo65GuObnXvFq0rKNsbJy5Co8kYSrsEZMd8jNrFDXQFG4Sp+r8EQSrsIaDW8Cux/sW3zFq5AjV+HJJFyFNXyDoMENsHfRFa8iNdvocw31D3VWVUI4jYSrsE5Md+OxLyf3XdHiDpsDgILCAmdWJYRTSLgK6zT7MygbbPniihYvGiVQ1D0ghCeRcBXWCa4FjXvArzOhIK/cixf1tRYNyRLCk0i4Cmu1ewDOJsOeBeVeNDzADNdsCVfheSRchbUa3QbBtWHTtHIvWjOoJgCJZxKdXJQQFSfhKqxls0Pb+2H/snI/FTbMP4yqvlVJOF2+5YRwBwlXYb22g42rtjZ/Vq7FlFI0CGnAgYwDLipMiCsn4SqsFxoNjW43T2zll2vRZmHN2Jm6U4ZjCY8j4So8Q7sHIPM47CvfRQWtI1tzNu8sBzMOuqYuIa6QhKvwDDHdoGrNcp/Yah3RGoCtJ533uG4hnEHCVXgGu49xYmvfEkg/UubF6gXXI9g3mK0pEq7Cs0i4Cs9x3f3G+69lP7GllKJNVBs2ndjkoqKEuDISrsJzhNY1xr1u/qxcJ7Y61+pMwukEGe8qPIqEq/AssQ/CmWOw49syL9KllvHYmNVHK/bYGCGcScJVeJbGPaF6S1jxWpmPXusF16N2ldqsOrrKxcUJUXYSrsKz2Gxw8/PGrQi3zi7TIkopbqxzI2uT1nI276yLCxSibCRchedpeifUvBZ+/meZ75bVvX53cgpyWHFkhWtrE6KMJFyF51EKbnkR0g/Bls/LtEjbqLZEBUaxMGGhi4sTomwkXIVniukGddrDin9Cbum/6tuUjW71urH66GpO5552Q4FCXJ6Eq/BMSkG3CcbIgV/eLtMiPRr0IK8wj2WHl7m4OCFKJ+EqPFfdTtD6XlgzGdJKv/NV64jW1KlSh3n757mhOCEuT8JVeLbbxoHNAYteKnVWpRS9Y3qz8fhGDp0+5IbihLg0CVfh2YJrwk3/gD3/g/ilpc5+d6O7sSs73+4r+0UIQrhCpQ3XU2dzOXU21+oyhDN0egTCGsLC50odmhUVGMUNdW7g+/3fk1dY/oceCuEslTZcE9PP0X7iTwz5ZANzNh4hPUuC1mv5+EGP1+HkXlj7bqmz94npw8lzJ1l+eLkbihOiZJU2XBtFVeGhGxuScPIsz3y9ldgJP3H/1PV8sf4wJzNzrC5PlFfj7tD0z7D8NTi577Kz3lD7BmpXqc3nu8o2RlYIV3B5uCqlEpRS25RSW5RScWZbmFJqiVJqn/lezWxXSqnJSql4pdRWpdR1xdYz1Jx/n1JqaGnbDXDYebZHU37+x838MOpPDL+hAYfTsnjh2210mPgT905Zy7TVB0nKOOe6nRfOdedb4AiA70dBYeElZ7Pb7AxsOpDNyZvZkbrDjQUK8TultXbtBpRKAGK11ieLtb0BpGmtX1dKPQdU01o/q5S6A3gMuAPoCPxba91RKRUGxAGxgAY2Ae201qcutd3Y2FgdFxd3QZvWml1JZ1i44zgLtyex90QmAG2iQ+nZsgY9W9akbnigE/deON2WL+C7v0HPN6Djw5ec7UzuGW776ja61u3Kqze86sYCRWWnlNqktY4tdT6LwnUPcLPWOkkpVRNYobVuopSaYn7+svh8RS+t9cNm+wXzlaSkcL3Y/pRMFm4/zoLtSWw/alzV07xmsBG0rWrQKKrqFe+3cBGt4fO+cGgtPLIGqtW/5Kyvb3id2Xtms7jPYiIDI91Xo6jUyhqu7uhz1cBipdQmpdRIs6261joJwHyPMttrA8Wf8ZFotl2qvUKuiazCo7c0Yv5jN/DLM7fw4h3N8HfYeGvJXm57eyW3vf0zby3ew/ajGbj6h5AoI6Xgz++AzQ7fjLzsbQkHNh1IQWEBs/bMcmOBQhh83LCNLlrrY0qpKGCJUmr3ZeZVJbTpy7RfuLAR3iMB6tatW64io8MCeejGhjx0Y0OOZ2SzeOdxFmw7znvL4/nPsnjqhgXSo2UNerSsQZs6odhsJZUk3CI0Gu58G74ZAb+8CTc/V+JsdYPrcmvdW/ly95c82OJBqvhWcXOh4mrm8iNXrfUx8z0Z+BboAJwwuwMw35PN2ROB6GKL1wGOXab94m19qLWO1VrHRkZe+a+BNUL8GXJ9fb4c2YmNL97G671b0SAiiE9XH6T3f9fQ+fVljJ23g3UHUikolCNaS7TuB63vM25LeHjdJWd7qNVDnMk9w+w9Zbs3rBDO4tI+V6VUEGDTWp8xPy8BxgNdgdRiJ7TCtNbPKKXuBEbx+wmtyVrrDuYJrU1A0eiBzRgntNIute2y9LmWV8a5PJbtPsGCbcf5eW8KOfmFhAf5cluz6tzevDp/ionA32F36jbFZWSfhik3GCMH/rYK/ENKnO3hJQ+zO203i/oswt/H381FisrGI05oKaUaYhytgtEF8YXWeqJSKhyYA9QFDgP9tNZpSikFvAv0ALKAB7XWRcO3hgEvmOuaqLX+9HLbdkW4Fnc2J58Ve1JYuOM4K3YncyYnnwCHnZsaR9KtRXVubRpFaKCvy7YvTEc2wifdoUUv6POx0Sd7kY3HNzJs0TCe7/A8A5sNtKBIUZl4RLhaydXhWlxufiHrDqSyeOdxluw8wYnTOdhtig71w+jWwjiqrVNNhni5zM//guUT4K53f388dzFaa4YuHErS2SR+7PUjDrvDgiJFZSHh6sZwLa6wULPtaAZLdp5g8c7j58fSNq8ZzO3Nq9OtRXWa1wxGlXCEJa5QYQHM7G30vQ5fbDwi5iIrE1fy6NJHGd95PL1iellQpKgsJFwtCteLJZw8ez5o4w6dQmuoHRpwPmg71A/Dx15pr0J2n7MnYcqNYHfAyBUQUO2CyVpr7vvffWTkZPDDPT/I0au4YhKuHhKuxZ3MzGHZrmQW7zzOL/tOkpNfSEiAg65No7i9eXVuaBxJFT93jI6rpI5shE97QqPb4L4vjCfJFvNL4i88svQRRncaTf8m/S0qUng7CVcPDNfisnLzWbn3JIt3HmfZ7mTSs/Jw2BWdGoZzS5MoujaLol54kNVlep/1H8KCf8Cto+HGpy+YpLVmyIIhHDt7jB97/4if3c+iIoU3k3D18HAtLr+gkLhDp1i+O5llu5PZl2z0014TGUTXZsbIg3b1quGQ7oPSaQ3fPATb5sKAL6FJzwsmb0jawPDFw3m2/bMMbj7YoiKFN5Nw9aJwvdjh1CyW7T7B0t3JrD+QRm5BIVX9fbipcSRdm0VxU+MowoJkmNcl5WYZ3QOp8TB8CVRvfsHkEYtGsC99Hwt6LyDQIaM4RPlIuHpxuBaXmZPPqn0nWbb7BMt2p3AyMwebguvqVuOWpkb3QZPqVWX0wcUyjsJHt4CPPzy0HILCz0/akryF+xfczxPXPcHwVsMtLFJ4IwnXShKuxRUWarYfy2DpLqP7YNvRDMAYfXBj4whujImkc6MIQgLkTDgAiXHw6R1Qpz3c/y34/H60/8hPj/Bbym8s6LOAYN9gC4sU3kbCtRKG68VOnM5m+e5klu9JZk18Kmdy8rHbFG2jQ7mpcSQ3No6kVe2Qq/smM1vnGH2w7R6EP086fwXXrtRd9J/fn+Eth/NEuycsLlJ4EwnXqyBci8srKOTXw+ms3JvCz3tTzh/VVgt0cEOMEbQ3No4gqupVeG39T2Nh1STo/ipc/+j55ud/eZ4lh5Ywv9d8agTVsK4+4VUkXK+ycL1YamYOq+JP8vOeFFbuS+FkpvGAxmY1g82j2gja1auGn89VcKOZwkL4aijs+gH6TYMW9wBwLPMYf/n2L/Ro0IOJf5pobY3Ca0i4XuXhWlxhoWZn0mlW7kvh5z0pbDp0ivxCjb/DRvv6YXS+JoIujcJpUSsEe2XtQsg7BzPuhmNbYMj3UO96AN6Oe5tpO6bx1V++oklYE4uLFN5AwlXC9ZIyc/JZuz+V1fEnWbP/5Pn7HwT7+9CpYThdGkXQ+ZpwGkVVqVyjELLSYOrtkJVqDNGKiCEjJ4M7vrmDlhEtmXL7FKsrFF5AwlXCtcySz2Szdn8qa+JTWXPgJEfSjCfiRlX1o/M14XS+JoLOjcIrx5290g4aAesIgOE/QdXqzNgxg3/F/Yspt02hc+3OVlcoPJyEq4TrFTuSlsWa/SdZHZ/Kmv2pnMzMAaBWiD/tG4TRvn4YHRqE0SiyineORDi6GabdCRGN4YH/kevjy13f3UUVRxVm/3k2dttV0A8trpiEq4SrU2it2Zecydr9qWxISGPDwTRSzhhhWy3QQbt6YXRoUI329cNoWTvEey7R3bMQZg2AhrfAgFksOLKUZ1Y+I7ckFKWScJVwdQmtNYdSs9iQkMbGg2lsTEgjITULgACHnbZ1Q2lf3zi6vTY6hKr+HnxBw+bPYN4oaNGLwt4fcf/CBzh29hjze80nyCE3zRElK2u4yv3tRLkopagfEUT9iCD6xxrPjEw+nc3GhFNsOJjKhoRTTF62D62N8foxUVVoEx1K27rVaBMdSuPqVT1nRMJ190N2Oix+CZt/CM92fIZBCwbz0daP5MICUWFy5CqcLuNcHluOpLPlcDq/HjnFliPppGflARDoa6d1nRDaRFejTXQILWuHUDs0wNpRCUvHwy9vwZ/+jxf8clmYsJDv7/6e6ODoC2abOHEiU6ZM4ciRIzRt2pRp06bRsWPHElfZtGlTJk2aRM+ePUucXpoVK1bw5ptv0rdvX+Li4nj33XevaD3C+eTIVVgmJMDBTY0juamx8Xjzoq6EX4+cYsvhdLYcSWfqqgPkFRg/2EMDHbSsZQRty9rBtKwVQt2wQPedLLt1NJw7Basm8febn+Enmw9vbXqLd2555/wsa9euZf78+bRv35569epx/fXXEx0dfclVHj9+nJSUFHdULzyUhKtwueJdCb3a1gEgO6+AXUmn2XHsNDuOZbDtaMYFgVvVz4dmtYJpVqMqTWsG06RGVZpUr0qQK57UoBTc8SacS6f6ijcY0XEA/zm8lPVJ6+lY0zgyTUpKIjQ0lHXr1rF8+XLuuusu3njjDZYuXcqQIUNITU0lMDCQYcOG0alTJ86cOcPo0aN54YUX6N+/P9988w1xcXHMnz+fBQsWcOLECUaMGMFrr71GQkICAQEBxMTEMG3aNNasWcPLL79MVlYWNWvWxM9PburtjaRbQHiM3PxC9p44w45jGWw/aoTunuNnOJtbcH6eumGBNK5ehXrhQdQLD6RuWCD1woOoUy2g4iMV8nNh1gCy9y/jnsYtCQyKZM6f5+Bj8yEzM5MWLVqQnp7OoEGD+Pnnn/n444+56667qFOnDqtXr+avf/0rTZs25YUXXiA0NJTJkyfTs2dPrr/+evLz84mLi+P++++nT58+TJkyhRo1ajBt2jRCQkJ4/PHHqVq1Ktu2bWP79u1MnjyZDz/8kIyMDNq2bSvdAh5EugWE1/H1sZldAyHc295oKyzUHE0/x+7jZ9iddJrdJ84QfyKTVfEnyc4rPL+sTUGt0AAzcIOoHx54/nO98MCyHfH6+EL/z/Cf2Zsnj27jqch0vvn8OfqPnkuVw4dp4efHrb16cSYykiNHjvDiiy/i5+fHI488QmBgIEOHDuW99967YJWRkZE0bNiQbdu2kZaWxp49e7j22mtJT0/n6NGj3HjjjSQmJpKdnU1ISAgOh4OuXbsSGhqKzWbj3nvvZe/evc78YxZuIuEqPJrNpogOCyQ6LJDbm1c/3661JuVMDofSskg4eZbDaVkcSs3iUFoWC7cncco8gVYkoorv+aNc4/338I2o4vv7CTXfQBg4h9s/u4d25xJ5Va9kYMoWauqDnMhuz45Zs1Dh4TgcDtauXUt4eHipJ+PuvfdeNmzYwLx58+jVqxc5OcY44aFDh5KUlMTDDz/M448/TkJCAh06dKhclxxfxSRchVdSShEV7E9UsD/t64f9Yfrp7DwOpxYF7lkOp2aRkHqW9QdS+W7LUYr3hgX52ok+H7hG+CbkTqFXwiA2t8ii3QMvsOG/1bFzN69qPwYFreGlwYOZOnUqmZmZvPfeewwcOJDPPvuM9u2NQ2673U5WljH+t3fv3owaNYrPPvuMjz76iC+//JLQ0FDmzp1LTEwMtWvXJi0tjXfeeQdfX19WrFhB7969KSws5KuvvuLaa691y5+pMH5o5+QXkpVbwNmcfLJyC8jMySczJ5+z5ntZSbiKSinY33G+i+Fi2XkFJJ46x+G0s0b4pmZxOC2L+ORMlu9OIbfA6G4IsT9F19MTWN5+JSE1Ujmb7cuDhYrHT2URtXQdsXcM5Oie3zi0Zzuh1aoREBjIwMEPkJGVR40aNXjttdf44IMPWLt2LW3btmXjxo089dRTdOzYkaCgIJ5++mlGjx7NwIEDcTgc9OvXDx8fH8aMGcOoUaPIyspi8ODBFBQU/GEfrnZ5BYVk5xWQnVfIudwCzubmk5Wbb4ZiAefy8o13c9rv8xSQlVNAVl4BWWZ4Fi1X9LnQSaeh5ISWEMUUFGqSz2TToFkOtsAcqtfcS40HHiEmL5fUlc+zicY0D91FSlQdzmRf+ihGKSPgQwKMV3CAD4G+PgT62s+/B/naCfD1IcjPToDDTpCfDwG+dvx8bPjabTjsNnx9zPfzn9UFbe4arqa1pqBQU6ghv7CQvHxNbkEhecVeufn6988FheQVaPLyL/peUEhu/u/BmJ1fcP5zTl6B+b1o+u/z5BRvyy+koJwJ6OdjO/9nH+Rn/rn72i/4+zj/7mcn0GEn0M/4XsXPhyp+PgSZ73XDg+SElhDlZbcpaoYEUMM3gEP7IWF/dXTgE2QOfIeXW43H97//5tvX/WFQd7LzCkjPyiPjXB7pWblknMs7/zpd7HPGuTxOZ+eTdvbc70dJOflk5RVQ0WMbH5vCx66wKYUCbEqBMt5V0TtGN4rxHYwWKNTafBmfCwo1WmOG6O/Tyhtk5eGwK/x97Pg57Pg7bPgXvfsYYRcWZDOm+fxxetHnoqAM8DV/QJk/qAKLhacVVwVKuApRgokTYeRIyMqCQ0seoGmXH/lPrb18N+xhuPkHAPwddmqE2KkRcmWPztFak51X+IdfW3PzjaO7oqPAnPwLj/qKjgR/n0ejzSDUZlCC8V70XZvbKz7dblMopbApsCvjs91mfLfZjMC2qwu/2xQ4zKNqh48NX7v6/fsFn234+lz03W7D4aPwtRsh6edjw8dbbvRzBSRchSjBoEHG+4svwuHDdnIWjiN95EA+rhXOMzN6wf3fQJ1SfzO8LKUUAeYRl6h8Ku+PDSEqaNAgSEgwHsF1YE1L+jbpyxcBNvZWCYUZ98Dh9VaXKDyYhKsQZfR428ep6hvMxAYt0FUiYWZvOLTG6rKEh5JwFaKMQv1DeeK6J9icup35tz4JVWvCzD5w8BerSxMeSMJViHLoFdOLVhGteGvHVM4Mmg0h0fB5PziwwurShIeRcBWiHGzKxosdXyQtO43/xs+FB/4HYQ3gi3shfqnV5QkPIuEqRDm1iGhBv8b9+HL3l+zJTYOh8yE8Br4cAPuWWF2e8BASrkJcgceve5xg32AmrJtAYWA1GDoPIpvArIHGww/FVU/CVYgrEOIXwpOxT7IlZQvf7PsGAsOMgK3eAmYPhh3fWV2isJiEqxBX6O5r7ia2eiyTNk0i9VwqBFSDId9D7XYw90HY8oXVJQoLSbgKcYWUUoy+fjRZ+Vm8Gfem0egfYly91eBG+O5vsOEja4sUlpFwFaICGoY0ZFjLYcw/MJ91SeuMRt8gGDAbmtwJPz4Nv7xtbZHCEhKuQlTQQ60eIrpqNBPWTSCnwHjKAA5/6D8dWvWDpeOMx3dX0tt7ipJJuApRQf4+/rzU6SUOnT7EJ9s++X2C3QG9psB1Q+GXt2Dhc8aNCsRVQcJVCCfoXKszPRv05KNtH5GQkfD7BJsd/vJvuH4UrP8A5j0GhfJkgauBhKsQTvJM+2fwt/szYd0ELnjCh1LQbQLc9BxsmQlzhxmP8RaVmoSrEE4SERDB36/7O+uPr2f+gfkXTlQKbnneCNmd3xljYfPOWVOocAuvCVelVA+l1B6lVLxS6jmr6xGiJP2a9KN1RGvejHuTjJyMP87Q+TH48yTYtxg+6w3ZJcwjKgWvCFellB14D+gJNAcGKKWaW1uVEH9kUzZevv5lMnIymLRpUskzxQ6DvlMhcSN8eidkJru3SOEWXhGuQAcgXmt9QGudC8wC7ra4JiFK1CSsCYObDebrfV/za/KvJc/Usg8MnAVp++GT7nAqwa01CtfzlnCtDRwp9j3RbLuAUmqkUipOKRWXkpLituKEuNgjbR6hRlANxq8dT15hXskzNbrNuFw2Kw2rHSrCAAAYAklEQVSmdocTO91bpHApbwnXkp6L+4cR2VrrD7XWsVrr2MjISDeUJUTJAh2BvNjxReLT45m2fdqlZ4zuAMMWGie8Pu0hz+WqRLwlXBOB6GLf6wDHLKpFiDK5Ofpmbq93Ox/89sGFY18vFtUMhi2CwHCYcbfcE7aS8JZw3QjEKKUaKKV8gfuAeRbXJESpnu/wPH52P8avG3/h2NeLVatnBGxEDHx5H2yb674ihUt4RbhqrfOBUcAiYBcwR2u9w9qqhChdZGAkT8Y+ycbjG/k2/tvLz1wlCh6YD9Gd4OsRckctL+cV4Qqgtf5Ra91Ya32N1nqi1fUIUVa9Y3rTrno73ox7k5PnTl5+Zv8QGPw1NLnDuKPWitflhi9eymvCVQhvZVM2xlw/hpz8HF7f8HrpCzj8of8MaDMIVrwGC56RG754IQlXIdygQUgDHr72YRYlLGLFkRWlL2D3gbveNW74suFD+HoY5Oe4vE7hPBKuQrjJgy0epFFoIyasm0BmbmbpC9hs0H2icT+CHd/CzD5yuawXkXAVwk0cdgfjOo8jOSuZyb9OLvuCnR+D3h/B4bXG5bJnjruuSOE0Eq5CuFHryNYMaDqAWbtnsSV5SzkW7A8D50DaAfj4dji5z3VFCqeQcBXCzR6/7nGqB1Vn3Npx5BVc4tLYkjTqagzVysuCqd0gMc51RYoKk3AVws2CHEG81PEl4tPj+WT7J6UvUFzt62D4YmPI1vS/wN5FrilSVJiEqxAWuCn6JnrU78GUrVM4mHGwfAuHX2MEbERj+HIA/DrTNUWKCpFwFcIiz3Z4lgCfAMatHUehLuc41qKruRrcCN8/CivflIsNPIyEqxAWiQiI4OnYp9l0YhNf7/u6/Cvwq2qc5GrVH5a9Yl5sIA8/9BQSrkJY6J5G99ChRgcmxU0iJesK7kHs42s8vrvzY8bFBnMfhLxs5xcqyk3CVQgLKaV4+fqXyS3MZeL6K7xlhs1mXGjQbSLs/N642OBcunMLFeUm4SqExeoF1+ORNo+w9PBSlhyqwL1cO4+CPlPhyHr49A44Lbc8tpKEqxAeYEjzITQPb87EdRNLfmpsWbXqC4O+gvRDxsUGybucV6QoFwlXITyAj82H8Z3Hk5GTwRsb36jYyq65BR78EQrzjWdzHVzpnCJFuUi4CuEhmoQ1YVirYczbP49VR1dVbGU1r4URS6BqDfisN2z9yjlFijKTcBXCgzzc+mEahjRk/NrxnM07W7GVhdaF4YsguiN8MwJ+eVvGwrqRhKsQHsTX7su4zuM4fvY472x6p+IrDKgG938DLfvA0nHwv6egIL/i6xWlknAVwsO0iWrDoGaDmLVnFptObKr4Cn38oPfH0OUJiJsKswdDbgWPikWpJFyF8ECPtX2M2lVqM2bNGLLznXBRgM0Gt4+DO96EfYtg2p8h8wouWhBlJuEqhAcKdAQytvNYDp0+xPu/ve+8FXd4CO793BiiNfU2OBnvvHWLC0i4CuGhOtXsRO+Y3kzfMZ0dqU58knzTO4ybvuRkwtTb4fB6561bnCfhKoQHeyr2KcL8wxizegx5heW4sXZp6sQaQ7UCQmHGXbDrB+etWwASrkJ4tGDfYF7q9BJ7Tu3h0+2fOnflYQ1h+BKo0Qpm3w/rPnDu+q9yEq5CeLhb695Kj/o9+OC3D9ifvt+5Kw+KgCHzoOmdsPBZWPgCFJbz3rKiRBKuQniB5zo8R5AjiDFrxlDg7Hu2+gZC/xnQ4WFY957cttBJJFyF8ALhAeE82+FZfkv5jS93f+n8Ddjs0POfxq0Ld34Hn90DWWnO385VRMJVCC9xZ4M7uaH2DUz+dTKJZxKdvwGljJtu9/0Ejm4ynjB7KsH527lKSLgK4SWKbqxtUzbGrR2HdtV9Alr2gSHfw9kU+Pg2eYT3FZJwFcKL1AiqwZPtnmRd0jq+i//OdRuq19kYSeAbBNPuhJ3zXLetSkrCVQgv07dxX2Krx/Kvjf8iOSvZdRuKbAwjlhpDteYMgTXvyl21ykHCVQgvY1M2xnYeS25hLhPWTXBd9wAYQ7WG/gDN/gKLX4Qfn5a7apWRhKsQXqhecD1GtRnF8iPLWZiw0LUbcwRAv+nQ+XHY+DHMGmhcOisuS8JVCC91f/P7aRXRilfXv0rquVTXbsxmg26vwJ1vQ/wS+LQnnE5y7Ta9nISrEF7KbrPzSpdXOJt3llfXv+qejbYfDgPnQNoB+LgrHN/unu16IQlXIbzYNaHX8Ldr/8biQ4sr9lju8oi5HYYtNE5ufdID4n9yz3a9jISrEF7uwZYP0jy8ORPWTeBU9in3bLRGKxjxE1SrD5/3hzgn31SmEpBwFcLL+dh8eKXLK5zOPc3rG15334ZDasOwBcajvOc/AUvGyE1fipFwFaISaFytMSNbj+THgz+y7PAy923YryoMmA2xw2D1O/D1MLnpi0nCVYhKYkSrETSp1oRX1r1CRk6G+zZs9zFGEdw+HnZ8a9x8+6yLRy94AQlXISoJh83BK11eIT07nTc2vuHejSsFXf5ujIdN+s0YSXCVP59LwlWISqRZeDOGtxrOvP3zWJm40v0FtLgHhs6HnDPGAxAPrXF/DR5CwlWISubh1g/TKLQR49aO43TuafcXEN3eGEkQGAEz7oatX7m/Bg8g4SpEJeOwO5jQZQKp51J5c+Ob1hQR1gCGL4Y67eGbEbDyX1fdTV8kXIWohFpEtOCBFg/wbfy3rD662poiAsPg/m+hVX9YNgHmjYICJz7B1sO5LFyVUmOVUkeVUlvM1x3Fpj2vlIpXSu1RSnUv1t7DbItXSj1XrL2BUmq9UmqfUmq2UsrXVXULUVn8rc3faBjSkLFrx5KZa9GNVnz8oPeHcNOz8OtMmNkHzqVbU4ubufrIdZLWuo35+hFAKdUcuA9oAfQA/quUsiul7MB7QE+gOTDAnBfgn+a6YoBTwHAX1y2E1/Oz+zG+y3iSs5J5e9Pb1hWiFNzyAtzzvnGCa+rtkHbQunrcxIpugbuBWVrrHK31QSAe6GC+4rXWB7TWucAs4G6llAJuBeaay08H7rGgbiG8zrWR1zKk+RC+2vsV65LWWVtMm4Ew5DvITDaGah1eb209LubqcB2llNqqlPpEKVXNbKsNHCk2T6LZdqn2cCBda51/UfsfKKVGKqXilFJxKSkpztwPIbzWo20epX5wfcauGUtWXpa1xdT/k/F0A/9QmP4X2Da39GW8VIXCVSn1k1Jqewmvu4H3gWuANkAS8FbRYiWsSl9B+x8btf5Qax2rtY6NjIws9/4IURn5+/gzvst4jmUeY9KmSVaXAxGNjKFadWLh6+Gw4p+VciSBT0UW1lrfVpb5lFIfAfPNr4lAdLHJdYBj5ueS2k8CoUopH/Potfj8QogyaBvVlkHNBjFz10y61e9G+xrtrS2oaCTBD0/AilchbT/c9R/jBFgl4crRAjWLfe0FFN1Vdx5wn1LKTynVAIgBNgAbgRhzZIAvxkmvedp4QNByoK+5/FDge1fVLURl9Vjbx4iuGs2YNWM4l3/O6nKMIL3nv3DrS7B1tnHBQSW6J4Er+1zfUEptU0ptBW4B/g9Aa70DmAPsBBYCj2qtC8yj0lHAImAXMMecF+BZ4EmlVDxGH+xUF9YtRKUU6AhkXOdxHDlzhMmbJ1tdjkEpuPEf0PdTOLrZvCfBPqurcgrl0idHWig2NlbHxcVZXYYQHmfiuonM3jOb6T2n0zaqrdXl/O7IRvjyPijMg3tnQoMbra6oREqpTVrr2NLmkyu0hLjK/F+7/6NWlVq8vPplsvM96N6r0e3hoaVQtSZ81su46MCLSbgKcZUJdAQytvNYEk4n8N6W96wu50LV6hv3JKh/A3z/KPw01mufbiDhKsRVqFPNTvRt3JcZO2ewNWWr1eVcyD8EBn0F7R6EVZPgq6GQa/H43Csg4SrEVeqpdk8RFRjF6NWjySnIsbqcC9kd8OdJ0G0i7PoBpt0JZ05YXVW5SLgKcZWq4luFsdeP5UDGAT747QOry/kjpaDzKLjvc0jZbYwkOLGj9OU8hISrEFexLrW70KtRLz7Z/gnbT24vfQErNL0THlwAhfkwtTvs+8nqispEwlWIq9zT7Z8mIiCCl1a95HndA0VqtTHuSRBWH77oBxs+srqiUkm4CnGVC/YNZlzncezP2M/7W963upxLC6kNDy6EmO7w49Ow4DkoLLC6qkuScBVC8Kfaf6J3TG8+3fEp21K2WV3OpflVMfpgOz0K69+HWQONhyF6IAlXIQQAT8c+TWRAJC+uftFzuwcAbHbo8Src+RbsWwKf9ISMo1ZX9QcSrkIIAKr6VmV85/EczDjoeRcXlKT9CBg0B04lwEe3wrFfra7oAhKuQojzOtfuTJ+YPkzfMZ3fUn6zupzSNbrNuKLL7oBP74Bd80tfxk0kXIUQF3g69mmqB1Zn9OrRnnXvgUup3twYSRDVDGYPhtWTPeLm2xKuQogLVPGtwtjOY72newCganV44H/Q/G5YMhp+eBzycy0tScJVCPEHnWt1pl/jfkzfMZ0tyVusLqdsHAHGfWFv/AdsngEze0NWmmXlSLgKIUr0VOxT1Ayq6T3dAwA2m/Fkg15T4Mh64zHeqfutKcWSrQohPF6QI4hxXcaRcDqB//z6H6vLKZ9r74Mh8+DcKeOeBAmr3F6ChKsQ4pI61exE/8b9+WznZ/ya7FlDnUpV73rjRFdQFMy4x+0335ZwFUJc1pOxT1KrSi1Grx7tGQ82LI+wBubNt/9k3Hx7yctuu/m2hKsQ4rKCHEGM7zyeQ6cPeV/3AEBAKAyaC7HDYfW/Yc79kHvW5ZuVcBVClKpDzQ7c2+ReZu6cyaYTm6wup/zsPsblsj3fgD0/wic94PQxl25SwlUIUSZPtnvy/IMNva57AIybb3d8GAbMhrSDLr9kVsJVCFEmgY5AXunyCofPHGby5slWl3PlGneD4YvAVnTJ7A8u2YyEqxCizNrXaM+ApgOYuWsmccfjrC7nylVvYTzGu3oL45LZVZOcfsmshKsQolyeuO4J6lSpw+jVo8nK876nsp5XJQqG/gAt+xiP8P5+lFMvmZVwFUKUS6AjkPFdxpOYmci/N//b6nIqxhEAfabCzc/Dlpnw2T1Ou2RWwlUIUW7ta7RnULNBfLH7CzYe32h1ORWjFNz8nBGyiXHw0xjnrFZ7wK25XCE2NlbHxXlxn5AQHi4rL4u+P/SlUBfyzV3fEOgItLqkiju6CcKuMcbGXoJSapPWOra0VcmRqxDiihSNHjiWeYxJmyZZXY5z1G532WAtDwlXIcQVa1e9HYOaDWLWnllsSNpgdTkeRcJVCFEhj1/3OPWC6/Hympe9e/SAk0m4CiEqJMAn4Hz3wNub3ra6HI8h4SqEqLC2UW0Z3Hwws/fMZn3SeqvL8QgSrkIIp3is7WPUD67Py6tf5mye6+865ekkXIUQTlHUPZB0Nom34t6yuhzLSbgKIZymTVQbhjQfwld7v2LtsbVWl2MpCVchhFONajuK+sH1GbNmDJm5mVaXYxkJVyGEU/n7+DPhTxM4kXWCtzZdvd0DEq5CCKe7NvJahjYfyty9c1lzdI3V5VhCwlUI4RKPtn2UBiENGLN2DGdyz1hdjttJuAohXMLP7seELhNIzkrmzbg3rS7H7SRchRAu0zqyNQ+0eIBv9n3DqqOrrC7HrSRchRAu9WibR2kU2ogxa8ZwOve01eW4jYSrEMKlfO2+TOgygdRzqfxzwz+tLsdtJFyFEC7XIqIFw1sNZ97+eaw4ssLqctyiQuGqlOqnlNqhlCpUSsVeNO15pVS8UmqPUqp7sfYeZlu8Uuq5Yu0NlFLrlVL7lFKzlVK+Zruf+T3enF6/IjULIazx19Z/pXG1xoxbO46MnAyry3G5ih65bgd6AyuLNyqlmgP3AS2AHsB/lVJ2pZQdeA/oCTQHBpjzAvwTmKS1jgFOAcPN9uHAKa11I2CSOZ8Qwss47A4mdJlAenY6r214zepyXK5C4aq13qW13lPCpLuBWVrrHK31QSAe6GC+4rXWB7TWucAs4G6llAJuBeaay08H7im2runm57lAV3N+IYSXaRbejJGtR/K/A/9j6aGlVpfjUq7qc60NHCn2PdFsu1R7OJCutc6/qP2CdZnTM8z5hRBeaETrEbSJbEN6TrrVpbiUT2kzKKV+AmqUMOlFrfX3l1qshDZNyWGuLzP/5db1x40qNRIYCVC3bt1LlCaEsJLD5mB6z+nYVOU+n15quGqtb7uC9SYC0cW+1wGOmZ9Laj8JhCqlfMyj0+LzF60rUSnlA4QAaZeo9UPgQzAerX0FdQsh3KCyByu4rltgHnCfeaa/ARADbAA2AjHmyABfjJNe87TWGlgO9DWXHwp8X2xdQ83PfYFl5vxCCOGxKjoUq5dSKhG4HvifUmoRgNZ6BzAH2AksBB7VWheYR6WjgEXALmCOOS/As8CTSql4jD7VqWb7VCDcbH8SOD98SwghPJWqrAeBsbGxOi4uzuoyhBCVjFJqk9Y6trT5Kn/HhxBCWEDCVQghXEDCVQghXEDCVQghXEDCVQghXEDCVQghXEDCVQghXEDCVQghXEDCVQghXEDCVQghXEDCVQghXEDCVQghXEDCVQghXEDCVQghXEDCVQghXEDCVQghXEDCVQghXEDCVQghXEDCVQghXEDCVQghXEDCVQghXEDCVQghXKDSPlpbKZUCHHLjJiOAk27cnrtV5v2rzPsGsn/OVk9rHVnaTJU2XN1NKRVXlmeZe6vKvH+Ved9A9s8q0i0ghBAuIOEqhBAuIOHqPB9aXYCLVeb9q8z7BrJ/lpA+VyGEcAE5chVCCBeQcBVCCBeQcL0EpVQ/pdQOpVShUir2omnPK6XilVJ7lFLdi7X3MNvilVLPFWtvoJRar5Tap5SarZTyNdv9zO/x5vT67tq/4pRSY5VSR5VSW8zXHcWmOWVfPdWl9sMbKKUSlFLbzL+zOLMtTCm1xPzzX6KUqma2K6XUZHM/tyqlriu2nqHm/PuUUkMt2pdPlFLJSqntxdqcti9KqXbmn1W8uaxy+U5preVVwgtoBjQBVgCxxdqbA78BfkADYD9gN1/7gYaArzlPc3OZOcB95ucPgL+Znx8BPjA/3wfMtmhfxwJPl9DutH31xNfl9sMbXkACEHFR2xvAc+bn54B/mp/vABYACugErDfbw4AD5ns183M1C/blRuA6YLsr9gXYAFxvLrMA6OnqfZIj10vQWu/SWu8pYdLdwCytdY7W+iAQD3QwX/Fa6wNa61xgFnC3+RPyVmCuufx04J5i65pufp4LdHXLT9Syc+a+eqIS98Pimiqq+L+pi/+tzdCGdUCoUqom0B1YorVO01qfApYAPdxdtNZ6JZB2UbNT9sWcFqy1XquNpJ2BG/5dSriWX23gSLHviWbbpdrDgXStdf5F7Resy5yeYc5vhVHmr1ifFP36hXP31RNdaj+8hQYWK6U2KaVGmm3VtdZJAOZ7lNle3r9LT+Csfaltfr643aV8XL0BT6aU+gmoUcKkF7XW319qsRLaNCX/oNKXmf9y63K6y+0r8D7wirntV4C3gGGXqe9K9tUTeVu9F+uitT6mlIoCliildl9m3kvtqzf+GZR3XyzZx6s6XLXWt13BYolAdLHvdYBj5ueS2k9i/NriYx7RFZ+/aF2JSikfIIQ//mrkFGXdV6XUR8D8i+orUpF99USX2z+Pp7U+Zr4nK6W+xejmOKGUqqm1TjJ/HU42Z7/UviYCN1/UvsLFpZeVs/Yl0fx88fwuJd0C5TcPuM88098AiMHoLN8IxJhny30xTlDNM/t4lgN9zeWHAt8XW1fRGc2+wDJzfrcy/+EW6QUUnbF15r56ohL3w+KaykQpFaSUqlr0GeiG8fdW/N/Uxf/Whphn2jsBGeav2ouAbkqpamZ3UDezzRM4ZV/MaWeUUp3M8wJDcMe/S3efFfSWF0bIJAI5wAnzL6lo2osYZ5n3UOysI8ZZzL3mtBeLtTfECKV44CvAz2z3N7/Hm9MbWrSvnwHbgK0Y/3BrOntfPfV1qf3w9Jf55/yb+dpRVDtGv/dSYJ/5Hma2K+A9cz+3ceEImGHm31c88KBF+/MlkATkmf/vhjtzX4BYjB8++4F3Ma9OdeVLLn8VQggXkG4BIYRwAQlXIYRwAQlXIYRwAQlXIYRwAQlXIYRwAQlXIYRwAQlXIYRwgf8Hb7sKm0L/vxYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Solve the location of the transmitter.\n",
    "\n",
    "c = np.argmin(rec_times)\n",
    "p_c = np.expand_dims(towers[c], axis=0)\n",
    "t_c = rec_times[c]\n",
    "\n",
    "# Remove the c tower to allow for vectorization.\n",
    "all_p_i = np.delete(towers, c, axis=0)\n",
    "all_t_i = np.delete(rec_times, c, axis=0)\n",
    "\n",
    "def eval_solution(x):\n",
    "    \"\"\" x is 2 element array of x, y of the transmitter\"\"\"\n",
    "    return (\n",
    "          np.linalg.norm(x - p_c, axis=1)\n",
    "        - np.linalg.norm(x - all_p_i, axis=1) \n",
    "        + v*(all_t_i - t_c) \n",
    "    )\n",
    "\n",
    "# Initial guess.\n",
    "x_init = [0, 0]\n",
    "\n",
    "# Find a value of x such that eval_solution is minimized.\n",
    "# Remember the receive times have error added to them: rec_time_noise_stdd.\n",
    "res = least_squares(eval_solution, x_init)\n",
    "\n",
    "print(f\"Actual emitter location:    ({tx[0]:.1f}, {tx[1]:.1f}) \")\n",
    "print(f\"Calculated emitter locaion: ({res.x[0]:.1f}, {res.x[1]:.1f})\")\n",
    "print(f\"Error in metres:            {np.linalg.norm(tx-res.x):.1f}\")\n",
    "\n",
    "# And now plot the solution.\n",
    "fig, ax = plt.subplots(figsize=(5,5))\n",
    "max_width = max(tx_square_side, rx_square_side)/2\n",
    "ax.set_ylim((max_width*-1, max_width))\n",
    "ax.set_xlim((max_width*-1, max_width))\n",
    "\n",
    "for locus in loci:\n",
    "    ax.plot(locus[0], locus[1])\n",
    "\n",
    "ax.scatter(tx[0], tx[1], color='red')\n",
    "ax.annotate('Actual', (tx[0], tx[1]))\n",
    "\n",
    "ax.scatter(res.x[0], res.x[1], color='blue')\n",
    "ax.annotate('Solved', (res.x[0], res.x[1]))\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Animations\n",
    "Below cells generate the animations used above. Some variable re-use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from matplotlib.animation import FuncAnimation\n",
    "from IPython.display import HTML, Image\n",
    "from multilaterate import get_locus\n",
    "\n",
    "# For animation GIF output:\n",
    "# brew install imagemagick\n",
    "# brew works on mac, different for different OS I guess.\n",
    "\n",
    "# If using to_html5_video:\n",
    "'''\n",
    "plt.rcParams['animation.ffmpeg_path'] = '/usr/local/Cellar/' \\\n",
    "                                        'ffmpeg/4.0.2/bin/ffmpeg'\n",
    "'''\n",
    "\n",
    "\n",
    "towers = np.array([[5e3, 13e3], [16e3, 5e3], [5e3, 7e3],\n",
    "                   [12.5e3, 1e3]])\n",
    "tx = np.array([15e3, 8e3])\n",
    "distances = np.array([ ( (x[0]-tx[0])**2 + (x[1]-tx[1])**2 )**0.5\n",
    "                       for x in towers])\n",
    "rec_times = distances/v\n",
    "\n",
    "v = 3e8\n",
    "delta_d = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Create a visualisation for TDOA\n",
    "\n",
    "# Plot towers and transmission location.\n",
    "fig, ax = plt.subplots(figsize=(5,5))\n",
    "ax.set_ylim((0, 20e3))\n",
    "ax.set_xlim((0, 20e3))\n",
    "for i in range(2):\n",
    "    x = towers[i][0]\n",
    "    y = towers[i][1]\n",
    "    ax.scatter(x, y)\n",
    "    ax.annotate('Tower '+str(i), (x, y))\n",
    "ax.scatter(tx[0], tx[1])\n",
    "ax.annotate('Tx', (tx[0], tx[1]))\n",
    "\n",
    "circle = plt.Circle((tx[0], tx[1]),\n",
    "                    radius=1, fill=False)\n",
    "ax.add_artist(circle)\n",
    "\n",
    "# Annotations, to be updated during animation\n",
    "cur_time = ax.annotate('t = 0', (1e3, 19e3))\n",
    "t0 = ax.annotate('Tower 0 received at t = ', (1e3, 18e3))\n",
    "t1 = ax.annotate('Tower 1 received at t = ', (1e3, 17e3))\n",
    "TDOA_ann = ax.annotate('TDOA = ', (1e3, 16e3))\n",
    "\n",
    "v_vec = ax.arrow(tx[0], tx[1], 0, 1e3,\n",
    "                 head_width=500, head_length=500, fc='k', ec='k')\n",
    "v_ann = ax.annotate('v = {:.0E} m/s'.format(v), (tx[0], tx[1]+1e3))\n",
    "\n",
    "n_frames = 300\n",
    "max_seconds = 50e-6\n",
    "\n",
    "t0_rec = 0\n",
    "t1_rec = 0\n",
    "TDOA = 0\n",
    "\n",
    "def animate(i):\n",
    "    global t0_rec, t1_rec, TDOA, v_vec, v_ann\n",
    "    \n",
    "    t = i/n_frames * max_seconds\n",
    "    radius = v * t\n",
    "    \n",
    "    v_vec.remove()\n",
    "    v_vec = ax.arrow(tx[0], tx[1]+radius, 0, 1e3,\n",
    "                 head_width=500, head_length=500, fc='k', ec='k')\n",
    "    v_ann.set_position((tx[0] - 0.5e3, tx[1]+radius+1.5e3))\n",
    "\n",
    "    circle.radius = radius\n",
    "    \n",
    "    cur_time.set_text('t = {:.2E} s'.format(t))\n",
    "    \n",
    "    if(t >= rec_times[0] and t0_rec == 0):\n",
    "        t0_rec = t\n",
    "        t0.set_text(t0.get_text() + '{:.2E} s'.format(t0_rec))\n",
    "    if(t >= rec_times[1] and t1_rec == 0):\n",
    "        t1_rec = t\n",
    "        t1.set_text(t1.get_text() + '{:.2E} s'.format(t1_rec))\n",
    "    if(t0_rec != 0 and t1_rec != 0 and TDOA == 0):\n",
    "        TDOA = abs(t1_rec-t0_rec)\n",
    "        TDOA_ann.set_text(TDOA_ann.get_text() \\\n",
    "                          + '{:.2E} s'.format(TDOA))\n",
    "    \n",
    "anim = FuncAnimation(fig, animate,\n",
    "                     frames=n_frames, interval=16.6, blit=False)\n",
    "\n",
    "# HTML(anim.to_html5_video()) # doesn't play nice with github markdown\n",
    "anim.save('animations/tdoa.gif', writer='imagemagick', fps=60)\n",
    "plt.close()\n",
    "Image(url='animations/tdoa.gif')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "## Create a visualisation for locus\n",
    "\n",
    "TDOA_dist = v*TDOA\n",
    "\n",
    "# Plot towers and transmission location.\n",
    "fig, ax = plt.subplots(figsize=(5,5))\n",
    "ax.set_ylim((0, 20e3))\n",
    "ax.set_xlim((0, 20e3))\n",
    "for i in range(2):\n",
    "    x = towers[i][0]\n",
    "    y = towers[i][1]\n",
    "    ax.scatter(x, y)\n",
    "    ax.annotate('Tower '+str(i), (x, y))\n",
    "ax.scatter(tx[0], tx[1])\n",
    "ax.annotate('Tx', (tx[0], tx[1]))\n",
    "\n",
    "circle0 = plt.Circle((towers[0][0], towers[0][1]),\n",
    "                    radius=1e3, fill=False)\n",
    "circle1 = plt.Circle((towers[1][0], towers[1][1]),\n",
    "                    radius=1e3, fill=False)\n",
    "ax.add_artist(circle0)\n",
    "ax.add_artist(circle1)\n",
    "\n",
    "# Annotations, to be updated during animation\n",
    "cur_d_ann = ax.annotate('d = 0', (1e3, 19e3))\n",
    "r0_ann = ax.annotate('Tower 0 radius r0 = ', (1e3, 18e3))\n",
    "r1_ann = ax.annotate('Tower 1 radius r1 = ', (1e3, 17e3))\n",
    "\n",
    "n_frames = 300\n",
    "max_d = 10e3\n",
    "\n",
    "locus_plot, = ax.plot([], [], marker=None)\n",
    "\n",
    "def animate(i):\n",
    "    # global t0_rec, t1_rec, TDOA, v_vec, v_ann\n",
    "    \n",
    "    d = i/n_frames * max_d\n",
    "\n",
    "    circle0.radius = d + TDOA_dist\n",
    "    circle1.radius = d\n",
    "    \n",
    "    cur_d_ann.set_text('d = {:.2E} m'.format(d))\n",
    "    r0_ann.set_text('Tower 0 radius r0 = {:.2E} m'.format(d \n",
    "                                                + TDOA_dist))\n",
    "    r1_ann.set_text('Tower 1 radius r1 = {:.2E} m'.format(d))\n",
    "    \n",
    "    locus = get_locus(tower_1 = (towers[0][0], towers[0][1]),\n",
    "                      tower_2 = (towers[1][0], towers[1][1]),\n",
    "                      time_1 = t0_rec, time_2 = t1_rec,\n",
    "                      v = v, delta_d = delta_d, max_d=d)\n",
    "    locus_plot.set_xdata(locus[0])\n",
    "    locus_plot.set_ydata(locus[1])\n",
    "\n",
    "anim = FuncAnimation(fig, animate,\n",
    "                     frames=n_frames, interval=16.6, blit=False)\n",
    "\n",
    "# HTML(anim.to_html5_video()) # doesn't play nice with github markdown\n",
    "anim.save('animations/locus.gif', writer='imagemagick', fps=60)\n",
    "plt.close()\n",
    "Image(url='animations/locus.gif')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Create a visualisation for loci\n",
    "\n",
    "# Plot towers and transmission location.\n",
    "fig, ax = plt.subplots(figsize=(5,5))\n",
    "ax.set_ylim((0, 20e3))\n",
    "ax.set_xlim((0, 20e3))\n",
    "\n",
    "for i in range(towers.shape[0]):\n",
    "    x = towers[i][0]\n",
    "    y = towers[i][1]\n",
    "    ax.scatter(x, y)\n",
    "    ax.annotate('Tower '+str(i), (x, y))\n",
    "ax.scatter(tx[0], tx[1])\n",
    "ax.annotate('Tx', (tx[0], tx[1]))\n",
    "\n",
    "circles = []\n",
    "locus_plots = []\n",
    "for i in range(towers.shape[0]):\n",
    "    circle = plt.Circle((towers[i][0], towers[i][1]),\n",
    "                         radius=0, fill=False)\n",
    "    ax.add_artist(circle)\n",
    "    circles.append(circle)\n",
    "    \n",
    "    locus_plot, = ax.plot([], [], marker=None)\n",
    "    locus_plots.append(locus_plot)\n",
    "\n",
    "# Annotations, to be updated during animation\n",
    "cur_d_ann = ax.annotate('d = 0', (1e3, 19e3))\n",
    "\n",
    "n_frames = 600\n",
    "max_d = 10e3\n",
    "\n",
    "def animate(i):\n",
    "    # global t0_rec, t1_rec, TDOA, v_vec, v_ann\n",
    "    \n",
    "    d = i/n_frames * max_d\n",
    "    cur_d_ann.set_text('d = {:.2E} m'.format(d))\n",
    "    \n",
    "    first_tower = int(np.argmin(rec_times))\n",
    "    circles[first_tower].radius = d\n",
    "    \n",
    "    for j in [x for x in range(towers.shape[0]) if x!= first_tower]:\n",
    "        # print('tower', str(first_tower), 'to', str(j))\n",
    "        locus = get_locus(tower_1=(towers[first_tower][0],\n",
    "                                   towers[first_tower][1]),\n",
    "                          tower_2=(towers[j][0], towers[j][1]),\n",
    "                          time_1=rec_times[first_tower],\n",
    "                          time_2=rec_times[j],\n",
    "                          v=v, delta_d=delta_d, max_d=d)\n",
    "        locus_plots[j].set_xdata(locus[0])\n",
    "        locus_plots[j].set_ydata(locus[1])\n",
    "        \n",
    "        TDOA_j = v * (rec_times[j] - rec_times[first_tower])\n",
    "        circles[j].radius = d + TDOA_j\n",
    "\n",
    "anim = FuncAnimation(fig, animate,\n",
    "                     frames=n_frames, interval=16.6, blit=False)\n",
    "\n",
    "# HTML(anim.to_html5_video()) # doesn't play nice with github markdown\n",
    "anim.save('animations/loci.gif', writer='imagemagick', fps=60)\n",
    "plt.close()\n",
    "Image(url='animations/loci.gif')"
   ]
  }
 ],
 "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}