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\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 = 3e8m/s. \n",
    "\n",
    "If this signal is received by two receivers (towers), the time difference of arrival (TDOA) between the first and second tower can be determined. TDOA = t<sub>1</sub> - t<sub>0</sub> where t<sub>0</sub> is the time the signal is received first and t<sub>1</sub> 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 TDOA to be calculated.\n",
    "\n",
    "![](animations/tdoa.gif)\n",
    "\n",
    "This 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\\*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).\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 r1=d, and the circle around tower 0 has a radius of r0=d+TDOA\\*v (noting TDOA\\*v is constant and TDOA comes from the previous animation). \n",
    "\n",
    "![](animations/locus.gif)\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",
    "![](animations/loci.gif)\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",
    "## 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"
   ]
  },
  {
   "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": "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
}