{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Survey Study: a tutorial\n",
    "\n",
    "This Jupyter notebook aims to help razorback users to compute impedance estimates from the data set shown in the paper in different ways for a two stage remote reference configuration:\n",
    "\n",
    "1- Ordinary Least Squares\n",
    "\n",
    "2- M-Estimator \n",
    "\n",
    "3- Bounded Influence \n",
    "\n",
    "This tutorial is designed for Metronix data format (.ats files).\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import razorback as rb # importing the razorback library\n",
    "import numpy as np # importing the numpy library as np\n",
    "import matplotlib.pyplot as plt # importing the matplotlib.pyplot library as plt\n",
    "%matplotlib inline\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# function for getting sensor information from .ats file header\n",
    "def sensor(ats_file):\n",
    "    header = rb.io.ats.read_ats_header(ats_file) #razorback function for ats file data importation\n",
    "    chan = header['channel_type'].decode()\n",
    "    stype = ''.join(c for c in header['sensor_type'].decode() if c.isprintable())\n",
    "    snum = header['sensor_serial_number']\n",
    "    sampling_rate = header['sampling_rate']\n",
    "    x1, y1, z1 = header['x1'], header['y1'], header['z1']\n",
    "    x2, y2, z2 = header['x2'], header['y2'], header['z2']\n",
    "    L = ((x1-x2)**2 + (y1-y2)**2 + (z1-z2)**2)**.5\n",
    "    return chan, L, stype, snum, sampling_rate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# function for getting calibration function from sensor information\n",
    "def calibration(ats_file, name_converter=None):\n",
    "    chan, L, stype, snum, sampling_rate = sensor(ats_file)\n",
    "    if chan in ('Ex', 'Ey'):\n",
    "        return L\n",
    "    elif chan in ('Hx', 'Hy', 'Hz'):\n",
    "        calib_name = f\"{stype}{snum:03d}.txt\"\n",
    "        if name_converter:\n",
    "            calib_name = name_converter.get(calib_name, calib_name)\n",
    "        return rb.calibrations.metronix(calib_name, sampling_rate)\n",
    "    raise Exception(f\"Unknown channel name: {chan}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The functions `sensor()` and `calibration()` will help to load the metadata of the data set (electric dipoles length, orientation, calibration files...).\n",
    "\n",
    "Other strategies for handling these metadata could be used, it's up to you to design your own.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import glob\n",
    "\n",
    "## Create the inventory containing all the time series\n",
    "\n",
    "# - files is the list of data files to load\n",
    "# - pattern is applied to each data file to extract the strings {site} and {channel}\n",
    "# - tag_template create the tag for the data file using {site} and {channel}\n",
    "files = glob.glob(\"data/site*/*/*.ats\")\n",
    "pattern = \"*/site{site}/*/*_T{channel}_*.ats\"\n",
    "tag_template = \"site{site}_{channel}\"\n",
    "\n",
    "# correcting incorrect information about calibration files in file headers\n",
    "name_converter = {\n",
    "    'UNKN_H104.txt': 'MFS07104.txt',\n",
    "    'UNKN_H105.txt': 'MFS07105.txt',\n",
    "}\n",
    "\n",
    "# creating and filling the inventory\n",
    "inv = rb.Inventory()\n",
    "for fname, [tag] in rb.utils.tags_from_path(files, pattern, tag_template):\n",
    "    calib = calibration(fname, name_converter)  # getting calibration for data file\n",
    "    signal = rb.io.ats.load_ats([fname], [calib])  # loading data file\n",
    "    inv.append(rb.SignalSet({tag:0}, signal))  # tagging and storing the signal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All the data are now loaded in `inventory`.\n",
    "We can use `inventory` to explore and handle the data set."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can check the number of files in your inventory `inv`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "24"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(inv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can display the tags/labels included in the inventory `inv`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'site002_Ex',\n",
       " 'site002_Ey',\n",
       " 'site002_Hx',\n",
       " 'site002_Hy',\n",
       " 'site002_Hz',\n",
       " 'site004_Ex',\n",
       " 'site004_Ey',\n",
       " 'site004_Hx',\n",
       " 'site004_Hy',\n",
       " 'site004_Hz',\n",
       " 'site006_Ex',\n",
       " 'site006_Ey',\n",
       " 'site006_Hx',\n",
       " 'site006_Hy',\n",
       " 'site006_Hz',\n",
       " 'site009_Ex',\n",
       " 'site009_Ey',\n",
       " 'site009_Hx',\n",
       " 'site009_Hy',\n",
       " 'site009_Hz',\n",
       " 'site099_Hx',\n",
       " 'site099_Hy',\n",
       " 'site100_Hx',\n",
       " 'site100_Hy'}"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inv.tags"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Creating (using `pack()`) and showing (using `print()`) the SignalSet object for site004 only:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SignalSet: 5 channels, 1 run\n",
      "tags: {'site004_Ex': (0,), 'site004_Ey': (1,), \n",
      "       'site004_Hx': (2,), 'site004_Hy': (3,), \n",
      "       'site004_Hz': (4,)}\n",
      "----------  -------------------  -------------------\n",
      "  sampling                start                 stop\n",
      "       128  2016-04-28 19:00:00  2016-04-29 07:00:00\n",
      "----------  -------------------  -------------------\n"
     ]
    }
   ],
   "source": [
    "print(inv.filter('site004*').pack())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Same operation for site099 (magnetic remote reference only):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SignalSet: 2 channels, 1 run\n",
      "tags: {'site099_Hx': (0,), 'site099_Hy': (1,)}\n",
      "----------  -------------------  -------------------\n",
      "  sampling                start                 stop\n",
      "       128  2016-04-28 19:00:00  2016-04-29 17:59:00\n",
      "----------  -------------------  -------------------\n"
     ]
    }
   ],
   "source": [
    "print(inv.filter('site099*').pack())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Creating and showing the SignalSet object content for the full inventory (including sites 002, 004, 006, 009, 100, 099).\n",
    "The full data set is reduced to maximal synchronous time section. The `pack()` function is narrowing the time range to the window of common synchronousness of the whole inventory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SignalSet: 24 channels, 1 run\n",
      "tags: {'site002_Ex': (0,), 'site002_Ey': (1,), \n",
      "       'site002_Hx': (2,), 'site002_Hy': (3,), \n",
      "       'site002_Hz': (4,), 'site004_Ex': (5,), \n",
      "       'site004_Ey': (6,), 'site004_Hx': (7,), \n",
      "       'site004_Hy': (8,), 'site004_Hz': (9,), \n",
      "       'site006_Ex': (10,), 'site006_Ey': (11,), \n",
      "       'site006_Hx': (12,), 'site006_Hy': (13,), \n",
      "       'site006_Hz': (14,), 'site009_Ex': (15,), \n",
      "       'site009_Ey': (16,), 'site009_Hx': (17,), \n",
      "       'site009_Hy': (18,), 'site009_Hz': (19,), \n",
      "       'site099_Hx': (20,), 'site099_Hy': (21,), \n",
      "       'site100_Hx': (22,), 'site100_Hy': (23,)}\n",
      "----------  -------------------  -------------------\n",
      "  sampling                start                 stop\n",
      "       128  2016-04-28 19:00:00  2016-04-29 04:00:05\n",
      "----------  -------------------  -------------------\n"
     ]
    }
   ],
   "source": [
    "print(inv.pack())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to prepare signal set from inventory to get it ready for TF estimation procedure\n",
    "from itertools import chain\n",
    "\n",
    "def prepare_signalset(inventory, local_site, remote_sites):\n",
    "    patterns = (f\"{e}*\" for e in [local_site, *remote_sites])\n",
    "    signalset = inventory.filter(*patterns).pack()\n",
    "    tags = signalset.tags\n",
    "    tags[\"E\"] = tags[f\"{local_site}_Ex\"] + tags[f\"{local_site}_Ey\"]\n",
    "    tags[\"B\"] = tags[f\"{local_site}_Hx\"] + tags[f\"{local_site}_Hy\"]\n",
    "    if remote_sites:\n",
    "        remote_names = tags.filter(*chain(*(\n",
    "            (f\"{e}_Hx\", f\"{e}_Hy\") for e in remote_sites\n",
    "        )))\n",
    "        tags[\"Bremote\"] = sum((tags[n] for n in remote_names), ())\n",
    "    return signalset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function `prepare_signalset()` build the SignalSet needed for processing a given `local_site` along with some `remote_sites`.\n",
    "\n",
    "The function starts by extracting the channels of interest and `pack` them in a SignalSet.\n",
    "Then that signaset is enriched with specific tags (`'E'`, `'B'` and `'Bremote'`) that will be used later in the TF estimate function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Showing the SignalSet object with `'E'`, `'B'` and `'Bremote'` tags for processing site004 using sites 100 and 099 as remote references"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SignalSet: 9 channels, 1 run\n",
      "tags: {'site004_Ex': (0,), 'site004_Ey': (1,), \n",
      "       'site004_Hx': (2,), 'site004_Hy': (3,), \n",
      "       'site004_Hz': (4,), 'site099_Hx': (5,), \n",
      "       'site099_Hy': (6,), 'site100_Hx': (7,), \n",
      "       'site100_Hy': (8,), 'B': (2, 3), \n",
      "       'E': (0, 1), 'Bremote': (8, 6, 5, 7)}\n",
      "----------  -------------------  -------------------\n",
      "  sampling                start                 stop\n",
      "       128  2016-04-28 19:00:00  2016-04-29 07:00:00\n",
      "----------  -------------------  -------------------\n"
     ]
    }
   ],
   "source": [
    "print(prepare_signalset(inv, 'site004', ['site100', 'site099']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Defining a frequency array in logscale for TF computation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.00000000e-03 1.39742149e-03 1.95278683e-03 2.72886629e-03\n",
      " 3.81337640e-03 5.32889414e-03 7.44671121e-03 1.04061943e-02\n",
      " 1.45418396e-02 2.03210792e-02 2.83971128e-02 3.96827358e-02\n",
      " 5.54535079e-02 7.74919238e-02 1.08288880e-01 1.51325208e-01\n",
      " 2.11465098e-01 2.95505873e-01 4.12946259e-01 5.77059977e-01\n",
      " 8.06396015e-01 1.12687512e+00 1.57471952e+00 2.20054690e+00\n",
      " 3.07509153e+00 4.29719900e+00 6.00499825e+00 8.39151362e+00\n",
      " 1.17264815e+01 1.63868373e+01 2.28993186e+01 3.20000000e+01]\n"
     ]
    }
   ],
   "source": [
    "# Definining your output frequency in logscale / you can reduce nb_freq if you want to make a quick test\n",
    "# as sampling frequency is 128, we go up to half a nyquist frequency which is 32 Hz\n",
    "# recordings are long enough to try to reach 1 mHz\n",
    "nb_freq=32\n",
    "freq = np.logspace(-3, np.log10(32), nb_freq)\n",
    "print(freq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computing two-stage OLS Impedance estimate for site004 with sites 100 and 99 as remote references "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First stage: a regression is performed to estimate the TF between magnetic field at site 4 and magnetic field at (sites 99 + 100) . Second Stage: the first stage TF is used to produce a synthetic magnetic field and a second regression is operated between the latter and site 4 electric field.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SignalSet: 9 channels, 1 run\n",
      "tags: {'site004_Ex': (0,), 'site004_Ey': (1,), \n",
      "       'site004_Hx': (2,), 'site004_Hy': (3,), \n",
      "       'site004_Hz': (4,), 'site099_Hx': (5,), \n",
      "       'site099_Hy': (6,), 'site100_Hx': (7,), \n",
      "       'site100_Hy': (8,), 'B': (2, 3), \n",
      "       'E': (0, 1), 'Bremote': (8, 6, 5, 7)}\n",
      "----------  -------------------  -------------------\n",
      "  sampling                start                 stop\n",
      "       128  2016-04-28 19:00:00  2016-04-29 07:00:00\n",
      "----------  -------------------  -------------------\n",
      "starting frequency 0.001\n",
      "starting frequency 0.00139742\n",
      "starting frequency 0.00195279\n",
      "starting frequency 0.00272887\n",
      "starting frequency 0.00381338\n",
      "starting frequency 0.00532889\n",
      "starting frequency 0.00744671\n",
      "starting frequency 0.0104062\n",
      "starting frequency 0.0145418\n",
      "starting frequency 0.0203211\n",
      "starting frequency 0.0283971\n",
      "starting frequency 0.0396827\n",
      "starting frequency 0.0554535\n",
      "starting frequency 0.0774919\n",
      "starting frequency 0.108289\n",
      "starting frequency 0.151325\n",
      "starting frequency 0.211465\n",
      "starting frequency 0.295506\n",
      "starting frequency 0.412946\n",
      "starting frequency 0.57706\n",
      "starting frequency 0.806396\n",
      "starting frequency 1.12688\n",
      "starting frequency 1.57472\n",
      "starting frequency 2.20055\n",
      "starting frequency 3.07509\n",
      "starting frequency 4.2972\n",
      "starting frequency 6.005\n",
      "starting frequency 8.39151\n",
      "starting frequency 11.7265\n",
      "starting frequency 16.3868\n",
      "starting frequency 22.8993\n",
      "starting frequency 32\n",
      "(32, 2, 2)\n"
     ]
    }
   ],
   "source": [
    "sig = prepare_signalset(inv, 'site004', ['site100', 'site099'])\n",
    "print(sig)\n",
    "ImpOLS = rb.utils.impedance(sig, freq ,remote='Bremote' )\n",
    "print(ImpOLS.impedance.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Showing typical apparent resistivity and phase results using matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-14-d5d0ebee0674>:5: RuntimeWarning: invalid value encountered in arcsin\n",
      "  rad_err = np.arcsin(res.error/abs(res.impedance))\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(-180.0, 180.0)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZsAAAEpCAYAAABFmo+GAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABgmklEQVR4nO2deZwcVbX4v2d6MguZMMAQ1qyssiQk7C2RtA6yKD6QKItgWA0gi7wfDyXKvkzwPZ8iCkiQGKJBXySIqKDgQLNNR1kyTGSHJCRDQhIiaWYCma3P749bnfRMeqnqrt5m7vfz6c9MVd2+99Ttqjp17j33HFFVLBaLxWLJJxXFFsBisVgsgx+rbCwWi8WSd6yysVgsFkvescrGYrFYLHnHKhuLxWKx5B2rbCwWi8WSd6yysfRDRM4UkceLLUepIiLjRERFpNKHuh4TkbP9kMtiKXWsshmCiMgUEWkRkaiI/FtEnheRwwBUdb6qHptQVkVkrxzaGiciT4nIJyLyhogcM+D4N0TkPRHZKCIPi8gOCcd+JCJvi0iH893p2cqRL0QkLCIXZPNdVT1BVe/Pos1qEbnP6bcOEVksIidkI4OI3CAiv8nmu7kiIueISJ+IdIrIxyLyioicmHA8rtg7nc9yEbk6Q50iIlc5182nIrJCRG4TkeqEMnNF5JYU3z9JRFodeT4UkWYRGefbSQ9hrLIZYojItsCfgZ8BOwC7AzcCXXlq8rfAYqAB+AHwoIiMdGQ5ALgH+CawM/AJcFfCdzcCXwHqgbOBn4rIZ9006oflUcJUAiuBqZi+uRZYUKYPxYiq1gHbYX7734nIdgPKbOeU+RpwrYh8MU19dwAzgOnACOAE4AvAgkyCOC9V84ArMf063pEp5uF8LKlQVfsZQh/gUGBDmuPnAM85/z8DKOah3wmc5uw/EWgFNgAtwMQUde2DUWIjEvY9C1zk/N8EPJBwbE+gO7H8gPoeAa5McSwEtAPfAz4Afo15mboaeBdYj3ng7OCUrwF+4+zfALwA7OwcWw4ck1D3DcBvnP/HOX1SCdwK9AGbnP75OSDAT4C1QBRoAw5MIXMYuCCx34EfAR8By4ATPPyubcC0NMe/B7wPdABvAo3A8U5/9zjyv+KUPRd43Sm7FLhwQF3fBVYDq4ALnP7YyzlW7ZzDCmAN8AugNtO15mxv49R12MC+TijzT+CqFPXt7fwehw/YP9q5Dr/gbM8Fbkny/a8BrcW+Rwfrx1o2Q4+3gD4RuV9EThCR7VMVVNWjnX8PUtU6Vf0/ETkYmANciLFW7gEeSRymSOAAYKmqdiTse8XZHz/+SkJ772IefvsMrEhEaoHDgFfTnNsuGGttLObt9nLgZIwFsBvmIX6nU/ZszNvraOc8LgI+TVP3VqjqDzDK81Knfy4FjgWOds5hO+A0jEJzwxEYRbAj8N/AfSIimb4kIjs77SXtGxHZF7gU8xAfARwHLFfVv2IU/v858h/kfGUt5oViW4zi+YnzuyMixwP/DzgG2AvTt4n80JFlknN8d+A6F+cQcNrqAd5LUeZI4EDgnRTVNALtqvrPxJ2quhJYBKSziABeBj4jIj8Rkc+LSF0muS3uscpmiKGqHwNTMG+M9wLrROQR54Hlhm8B96jqP1S1T82cQxdwZJKydZi3+0SimOENN8cT+QVGMf0tjWwx4HpV7VLVTzEK8Qeq2q6qXRgL5WvOEFsPRsns5ZzHS07f5EqPI/9nAFHV11V1tcvvvqeq96pqH3A/sCtmeDElIjIMmA/cr6pvpCjWh7E49heRYaq63FHsSVHVv6jqu2p4Gngc+Jxz+FTgV6r6qqp+ghmCjcsimOvjP1X1385LRhNweppTOFJENmCswx8BZ6nq2gFlPhSRT4EIZljr4RR17YixuJKx2jmeElVdirGQd8dYwR868ztW6fiAVTZDEOcBeI6qjsK8Ke4G3O7y62OBK0VkQ/yDsQ52S1K2E/N2nMi2mOEZN8cBEJH/ceQ8VVXTRY5dp6qbBsj6hwQ5X8c8eHfGDLP9DTNHsEpE/tt5cOeEqj6JGU67E1gjIrOdeTI3fJBQzyfOvykfdCJSgTmPbozlEt//WMKk+pmq+g5wBUbZrhWR34lIst8r/v0TRGSR4zyyAfgSWx7Uu2Hmi+Ik/j8SMxT2UkKf/9XZn4pFqrodsD1mmPRzScrsiOmH/8Iog1S/04cYBZ2MXZ3jaVHVRap6qqqOdGQ5GjPXaMkRq2yGOM7b8FzMw9wNK4FbVXW7hM82qvrbJGVfBfYQkURL5SC2DPe86mwDICJ7YN7A30rYdyNmkvdYF5bHQEW0EjPvkShrjaq+r6o9qnqjqu4PfBYzbBT3dtuIeWjG2cVDm6jqHap6CGaYcB/gqgxye8axIu7DKM5pqtqT0P4JzrBYnarOd/Y9oKpTMApYMcNdW8nvDIcuxFgZOzuK4FHMXBQYC2FUwldGJ/z/IWYo8oCE/q5XM7mfFlXtBL4NfFNEJic53qeq/4uxgL6dopongdEicviAcxqNsbybM8kxoM0XgIdwf29Y0mCVzRBDRD4jIleKyChnezRwBmZMOxlrgD0Stu8FLhKRIxw30+Ei8uUBCgUAVX0L40hwvYjUiMhXgYmYhxmY4Z+viMjnRGQ4cBPwUHyOR0RmAt8Avqiqbuc9EvkFcKuIjHXqGykiJzn/f15EJjhzBR9jhr/6nO+1AqeLyDARORQzcZyKfv0jIoc5fTMMo7Q2JdTrJ3cD+wFfcYYMUyIi+4rIFxxFsgmjEOIyrQHGOVYSQBVG4a8DesW4VB+bUN0C4FwR2U9EtiFhPkZVY5jr4ycispPT9u4icpybE3J+41+Sfo7nNuC7IlKT5PtvYX7z+SJypIgEHI/HhcDfVfXvCcUDzjUZ/1SJWRLwrQTZPwP8B6nvDYsXCuGFYD+l82HLePT7mIfh+5hJ/m2d4+fQ30PoIszb7AbMMBYYL6YXnH2rgd+T2oNsHMbr6lPM5PcxA45/A+O5tBH4I463mHNMMfNBnQmf76doJ4SZHE7cV4GZzH4TMzT3LtDkHDvD2b8R88C9A8frCaM8/uG09xfn2FbeaM52EGOJfeSUa8R4hnVi3vTnA3UpZA4zwBttwPHNXl4D9setk7gXXPxzZop2JmK8uDqAf2Nc33dzjjVgvOA+Al529l3i9MkGzDDd70jw3gJmYob8VgEXO7KMdo7VYOZplmKU+OvA5SnkSnbOo5zffOLAvnaOC8YivixFnRUYz7t3MNfcSoyzRU1CmblOvYmf5zAWzJ+cc+/EeCX+EBhW7Pt2MHzE6fyscN5GN6mZ0LRYLEMMEdkP+BdQraq9xZbHUrp4GkYTkQoxK77/IiJrgTeA1SLyqoj8j4jsnR8xLRZLqSAiX3WGnbbHvPn/ySoaSya8ztk8hVl4NxPYRVVHq+pOGK+NRcBtInKWzzJaLJbS4kLMnM67mLmfi4srjqUc8DSM5vjo9+RaxmKxWCxDi6znbKxSsVgsFotbsnJ9FpFfYhasrRSRf4jIvSJymc+ylQUisrOIPCMm+u7/isj3nf4pCUTkFyJybam2Lz5FHS72eRYbMRGUnyu2HJbsEJFbxESZ/iBz6fIk23U2n8Ms+BoNnAL8ARjum1TlxQyMi+u2qnqlqjapalYh53Ml2QNHVS9S1ZuLIc/A9kUkJCLt5dyOs77kSTHpGd5x1g4lHr/A2d8pIn9Nt1K/kIiJZPCmiMRE5Jwkx/9TRD5wzmuO9A/Jv4OI/EFMGoj3ROQbBRW+wPj1AuShvdGYSNP7q2q6BcRlTbbKZhEmvARqVmM/qqq3+SdWWTEWeE1z8SF3gQzukPllgfMb/BGzTmUHzIvGb0RkH+f4VMwak5Oc48swKRZKgVcwK+9fHnjAWXR5NWaN0DjMOqMbE4rciQmJszNwJnC3s1iyYDiLb8sSF/fuWGC9bh0Tzk3dkrAgt7TJZnEOcBRmsdZ/Yayc+mIvGHIp9xRMSPwNmMVe5zj76zF5LNZhIs5eA1Q4x84hReh3zOKwHsyN2ImJhHsDzgJAp8x0p871mLwjy3EWNjIg1DkDFiY6Zb+HWSTYhQlrHw+Z3wG8BnzVKbsfW1ard+KkEUjSxrcwC97+jYlFtVvCMcUs4nybLRGSJUk/1mAWzO3obF8D9LJlYegtwO2J7WMs308xwTLjixB3c/prgdP/HZgFe4em+P1Shu930U5W6QYGtH+gU58k7HscuNn5/0fAnQnHdnP6dM8U55MylD9bUiZc6ZzvauDchOMNzu/3MWbB5s0MWCCZos3ncK77hH0P4Cx2dbYbgQ+c/4fjROJOOP5r4LYU9e+JCRuzni2LWrdLOL5VqoMU9czFREl4FLPw9hinPxdi7tNlJCwWda6j3zu/YwewBBMqaKbTfysxIY8Sf5tHMPfBO8C3nP2p0i7UY8IDrXbkvwUIpJD9BuBBR5aPMWkYkn7fOa/E63WuU8eRbHlWvQKEEuoPY9JbPO98dy9M4NcnnPN5E2cBdkJf3olZoNyBWbC8Z8LxAxK+uwZn4TQ+3DP9+iXLh/YyzJvP1Zg3t1eBd7Opq1AfYIzT0WdgAvk1AJOcY/Mwb6wjMG92bwHnO8fOcS68bzkXx8WYldNx54q59H+Y38CW1eb7OxfQFEwYkB85dXlRNq2Y+FO1zr6vs+XheRrmRtw1QdaBK7I3t4FJIvUhcDAmJMnPgGcSyirmrX07p7/WAcen6M9ncPKnYB6477JFCT/DFiWY2H6/80vor02YYI8BYBYmOGOyNo8DXnLkE4yC3dVlO1dgLPJRzrnfA/zWOXYhZuX4No4Mh+AozgF1TGBrZfME8Afn//8F7ko4trvTpyelOJ8vYx7OggnV/wlwcMI59GJC+Axz+ucTYHvn+O8wN/9wjBJ8f+Bvn6LNZMrmFZxcRc72jo7cDcBk4NMB5f8Ls7YmWf17YUL5V2MCcD7DlhePfTEP/Xj0gnGkVsRzMS8UR2Gu9W2c3/46zL20B0ZBHzfgOjoO81I2D/Oc+oHTf98CliXU/zQmgnQNJh3COhzFx4AXRmffw841MxzYCaPgL0wh+w2Y+/xkR/badN9n6/t+d8xD/EvO97/obI90jocxUTcOcM613unXc53tgzH3+QEJfflv4HDn+Hzgd86xERgFeKXTFyOAI/y6Z/r1i9eHttPQs0n2VWdTV6E+mDecPyTZH8BYDfsn7LsQCDv/nwO8k3AsnuBpl4QfMpWyuS7+4yR8txtvyua8DOfVivMwI7OyuQ/474RjdZibYpyzrcCUhOMLgKtTtHszTogXTOiS72DiVg20ehLb73d+Cf3194Tt/RnwcEs49gXMi8CROJZnivNM1s7rJLxFY6IA9zjyn0eaJHAJ3xmGecB91/n/WOf3/JtzvBFzk0/EPGDuwbyxnuHyGn0Y+E7COXxK/1Ata51zDziyfybhWNPA3z5FG8mUzbskvFQ456YYZfA5HCsn4fi3cO4PF+2dDCx2/t/LOYdjyBACxvk95yVsHwGsSHJP/yrhOnoi4dhXMC8GAWd7hHNO22Fe3vron9RvFlusihvoPzqxM+YZUZuw7wzgqRSy30D/l7i03x94vWKsv18PqPNvwNnO/2HgpoRjpzHgmexce9cn9OUvE459CXgjQY7FKc4j53sm8ZPtWF+riHwncYeafCGlzGjMTTWQHTFvSokJm97DvF3E8RT6PYF+4did73oNKJkYwh0RmS4mR3o8hPuBZMjTMUCezeepJtLuelKcK+ZNOtV5Po25SQ7GDFk8gXk7PxKjnDOGc0/TZk2ycW7NLXz/WHJMN6DG1f9kjEXyAeZtcAFmuAtVbQauxwz1vId5WeiIHx+IpA/lD2YcP3Flfvz3GMmW1NBxkiYcc8nAVA/x/zuSHIsf7yAJIrKTmBQG74vIx5ihlh0B1GOqA/qf31hgN+mf2uL79M/3sybh/0+BD3VLKK14sNI6zH0Qz7cTZ+A9n8hYjAJendD2PRgLxa3sXr4/Fvj6gHOdQv/0CQPrP2JA+TPpH6081X2d6rkYr9e3FB3ZKpudMZF/V4nIn0XkVhH5epZ1FYqVmCGLgXyI0dZjE/aNwQxL5Eq/cOxisk02JBx3E8peE74/FhNV91KgQU3493+xJfy7bvXt/qwi4TzFxLZrILtzbcEMi3wVeFpVX8P025cxiigZmeTLiLoL35+snWzTDQxsv01Vp6pqg6oehxnO+WfC8TtVdW81kTUWYpTCvwbWI5lD+adjHWaILTG8/xgX30tFv1QPzv9r1ERhfguolP6hqA4idcbUWZj+n6iq2wJnkXBOmjrVQTISf8eVmGGwxN9vhKp+yd0p9mMVsIP0j1SeeM8PvH5WYiyTHRPa3lZV0zlJDJTdy/dXYiybxHMdrv2dsAbW//SA8nWq6iayQ6rnYvxYzvdMnKyUjZrkQvsB4zFDRW9hxgNLmfnAMSJyqohUikiDiExy3nwWYELRj3Ae6P8P80aWKw9iQuh/VkSqMPNciQ+TVuBLjmvpLpi3vnQMx1xk6wBE5Fz659pYA4xy2krGA5jw8JOch10T8A9VXe7prNhspb2EiRAcVy4tmCHIVMpmDdAgIvVe2wNP4fuTtZNtuoGBMkwUE5J+GxH5L8zb5lznWI2IHOh4CI0BZgM/VdWPklSVKZR/Spxr9iHgBkeO/TFprlMiJpZZDeb6G+bIGr//5wHni8j+YuKdXRM/J1Xd6LR1k5h0EkdhvO1+naKpETgOKiKyOwkvA5I+1UEm/gl8LCLfE5FaMekDDhSRw1x+fzNq0kS3ALOcfpgInI95RsCAtAtqMq0+DvyviGwrJkbkno73oZv2vH7/N5jnxnHOedaIcecflaL8n4F9ROSbYtJiDHPulf1ciPdnYBcRuUJEqp1n4BHOMV/umTg5ucypSb/7sqrer6q+J4jyE1VdgRmmuBIzWdbKlre5yzAPr6WYMe0HgDk+tPmqU/fvMFZOB2bMOj7k+GvM5OxyzMX4fxnqew0zCR3B3BATMB4pcZ7EvHF+ICJbDWM5wzzXYt6oV2PeaNKl7M3E05jhgX8mbI/ATAonk/8NjEPJUsc097oGZVuMZfcRWzz8fuSynZ9ivI8eF5EOzMRn/KbaBfNiEA+J/zSpXza+iem7tZg5mi8mDCHXYK6dTkyfRDD9vRXOEM7lmBedjzCpFh5x0wkOl2KGQj7AKIZfZSj/OObh/lmMEvwUk4USVf0rJgz/U5h+fQ8zHBjn25g5qLWYfr3YubaTcSNmaDWK8X56KOFYNWZe70NH7p0wQ2EZcRTsVzCT+cucOn6JmRzPhjMwc1KrMOsEr1fVJ5xjv3f+rheRuKv4dMwLwmuY3+tBUmcFTYbr7zvK8CRM36zDWBhXkeJ57VxLx2Lu5VWYvv0hpr/T4nz3i5i+/QDjifp557Bf9wyQQ7gai3fE5DLfAOytqsuKLI7FYrEUDN8WAznDQJYBiMhXnKGO4Zi38CUYS8ZisViGDH6uPL3Px7pcISZ0yC9E5EERKdUw5ydhTNtVwN7A6WrNSYvFMsTIJepzZaJbpoiMVdVc3C/j9czBeDasVdUDE/YfjxlDDGB8xm9LOFYB3Kuq5+favsVisVj8x7NlIyLfEpE3gZXO5OuTInIkZpLND+ZiQkYkthnArK04AbPo7wzHAwcR+Q/MpH6zT+1bLBaLxWe8poX+HsbrYaqq7uqsDbgN4yGUylfbE6r6DMZbLJHDMQsFl6pqN8a76ySn/COq+lnMIiaLxWKxlCBeIwmfC0zQhKRpqvq4iByDiTOUL3an/4rZdsyK2RAmxUE1ZkFcUkRkBiZCL8OHDz/kM5/5TN4EtVgslsHISy+99KGqjsz2+57D1muS7JyqukZEfpatEC5ItqpaVTWMiROUFlWdjVlbwKGHHqovvviir8JZLBbLYEdEcpqT9zpn866IfDmJEDeR3zmTdvqH5hiF8e6yWEqGSCTCrFmziEQixRbFYik5vFo23wYWisn09wpmBfOXMKvx3/RVsv68AOwtIuMx8YtOx6y4tlhKgkgkQmNjI93d3VRVVdHc3EwwGCy2WBZLyeDJsnFcmw/DrKnZiLEuzlLV6cD9fggkIr/FhPnYV0TaReR8x8X6UkyU0deBBWnCZVgsBSccDtPd3U1fXx/d3d2Ew+Fii2QpUSIRmDXL/B1KeLJsREScBYl/dT6bUdUfDiiTFap6Ror9j5LGCcBiKSahUIiqqqrNlk0oFCq2SJYSJBKBxkbo7oaqKmhuhqFiAHuds3lKRC5zItpuxoko+wURuZ8M0WctlsFIMBikubmZm2++2Q6hWVISDhtF09dn/g4lA9jrnM3xmAxtv3XmTzZgIt0GMFFlf6KqrX4KaLGUC8Fg0CoZS1pCIWPRxC2bXA3gSMQorFCo9C0kT8pGVTdh1tPc5eQU2RGTwndDHmSzWCyWQUUwaIbO/FAQ5TYk53mdTRxnvc1qH2WxWCwJRCIRwuEwoVDIWkyDiGDQH6WQbEiulC+TrJWNxWLJH9aVevASWRkhvDxMaFyI4Ojsf1O/h+TyjVU2FksJksyV2iqb8ieyMkLjvEa6+7qpClTRPL05a4XjZUiuFOZ2slY2Tlj/isQ0AxaLxR+sK/XgJLw8THdfN33aR3dfN+Hl4ZysGzdDcqUyt5NV8jQRuRRYA7wnIq+IyAX+imWxDG2sK/XgJDQuRFWgioAEqApUERoXynubpeJuna1lcyUm+vMHIrIr0CQio1T1Bv9Es1iGNtaVevARHB2keXqzL3M2bimVuZ2sMnWKyBLgIFWNOdsVQKuqTvRZPt+xUZ8tFks54Oc8ix91ichLqnpotjJka9ncDfxeRL6nqu8AY4BPshWiLCiFGTaLxTIk8HuexS9361zIStmo6l0i8gHwSxGZCGwLzBeRr2MsnLf9FLLolMoMm8ViGRKU2xoaN2TlIACgqg+pagjYCTgYeBL4LHCPP6KVEKUyw2axWIYE8XmWQCDDPIvbENIlEGo653U2jutzm/PxJc1AyeFlhs0Ot1ksQxO3976LcsEgNN++hPDC9YSmNRAMTkhej5sRlxIZmcnKQaCcydpBwM2F5OVHdVufVVxlQ7F+riWzZ7N+4UIapk1jwowZhWt4qJDhh+3p6aF9yRI2rVmDioAIsvPOUF29dV1dXbBmDaiCCORSLhqFDRu2bG+3HdTXb12X23IONTU1jBo1imHDhvXbXywHgaGHmxk2twOtbpRSibyNWNyRl5/LhfZaMns2e154IfsB3Y8/zhKwCgf8e5lz8cO2t7dTM7yGhpFj2ci21NFB1Y7bUD1q7Nb1rV4NvQnr4HfeGXbdNatyn3y0juql7yEKKtC1x1i22X7k1nV1dsJbb0EsBhUVsM8+UFeX9HRVlfXr19Pe3s748eOT90mWWGXjJ26H29woJS8zhNYCKjq+T+i61F7rFy5kP8yNrM42Q13Z+Pky5+KH3bRpE8OGb887HTsRo4IKYuys69g9mWwjRpgHfvzBP2JE8nNwUS5a2Uu0Aeq6oLMa6it72SZZXXV1RsF0dJh6UigaABGhoaGBdevWpSyTLVbZ+EkwyJLbb98ypJHqaeNGKblVXH4P3VmywveFcy61V8O0aXQ//jgK9Djbgxo317CfL3Muf9i+ih2IUQEIMYRe2T65bM6DP3TssRAIEH722bTl0imIEdUjWF1dQWdVjAqpYFR1CsUVry+NkklERFyV84pVNj4SiURovOIKE8/q2WdpnjAh+QpwNxH03EbZ83PozpI1fuYpAVw/5CbMmMESKP85Gz/nRP18mXP5wzZsV82HaxRVNdbBdknmYeLU1Zk2M5FBQdRV1bFPwz50dHUwonoEdVXulEmxyErZOLHR5qvqRz7LU9aEw2E+/fRTgMyRet3MAbkp4+fQXRwfvWoGDS49iHzrBg/aa8KMGe6GzoqxJN1PJeL2GvbzZS5eNkN/1dXBvvuKm5EqX6mrqit5JbMZVfX8AW4B3gEWYFJFSzb1FONzyCGHaL5oaWnR2tpaDQQCWltbqy0tLXlra0DDqk1N5m+6MrW1qoGA+ZuqbD7KZZKt1HF7rqVMS4tqRYUq5H4Ofl8jTU2mDJi/TU251VdgXnvtNc/fmTRpko4fP75wzwiPJDsn4EXN4dmb1aJOVb0G2Bu4DzgHeFtEmkRkTz8UYLniJVJvJBJh1qxZRPxYZBUMwsyZmd/Ompvh5pvTD6G5XcDqplz8jfXaa81fP861GIvTwmH49NOCL+r19VTDYeNCC7mfg5/XCLhfwej2Gi5xIpEIbW1tLFu2jMbGRn+eAQ533XUXBx54IGPHjuVnP/vZ5v2d3Z2s7lhNZ3enb215JhdNBRwE3A68gYmXthj471zqzPcnn5aNW1paWrSiokKBtBZQS0uLNjU1Ffbtx8+3VrdvrH7LFi/rl0VVhDdq35v00xL127Jx226RaFnRok3PNGnLiuSyebVsmpqaFOM8qIFAQJtyvS8cHnzwQT355JO1u7tbV61apQ0NDdrT06MdXR360qqX9IX3X9CXVr2kHV0dGevKh2WTrZK5HHgJ+BvwdWCYs78CeDcXgfL9KQVl09TUpIFAIO3FVrQhOdO4uxs/Uzm/HzbFHG4p8MPQbz2tqoX/vbyUK1FaVrRo7S21GrgxoLW31CZVOF6VjduXTa8cffTR+uabb27e3mWXXbSrq0tXfbxKX3j/hc2fVR+vylhXKSmbm4CxKY7tl4tAHuXYAzOU96Db75SCsnGjSNwopMT6Cm4BuaUYb8pentTFeBi6aLMo0xN50XDlTdMzTRq4MaDcgAZuDGjTM1v3SSnM2XR3d+t22223eXvVqlV64IEH6tq1a/XM6Wfqoy8+ql857Sv69yV/1zOnn6krV67Uc889V7u7u5PWV0rK5odu9mVZ9xxgLfCvAfuPB97EOCZcPeBYWSkb1cwKwq1lkw8LqODKy2/lkIdhHt8eDi0t2lI1VZvk+9pSNTWjwimoHizRCfhikg/LRlV16tSpOnXqVB8kNLS2tqqI6Lvvvqt9fX169tln6z333KOqqtddd52e8rVT9O1Vb2tHV4ded911euqpp2pHR+rhtFJSNi8n2deWiyAJ9RyNiSL9r4R9AeBdx5KpAl4B9k84nndlUwzrwU2bZvz3SIWrtaLiqJwtIL9NfFcPzGINezU1mVsgg5Lzs09aLrpfa9moAXq0lo3actH9WdeVF8p82Csf+D1no+q/svnVr36lZ511lh5++OF6wAEH6LXXXquqqh0dHXrWWWfpN77xjaTbqciHsvG0zkZELga+DewhIm0Jh0YAz3upKxWq+oyIjBuw+3DgHVVd6sjxO+Ak4DU/2sxEJBJhypQpxGIxamtrC5YT3k1a4IaGE4HvAFXEYt00NLybtFwkEqGxsdEsOK2qSnkO4XB48wriTGuFMi2hcL2ONBgkcvs/0ke4ddnm5nIECRMkBKQsFgoRqZpKuOcoQoHnCabwggqHw8RihwMhurqeTd0nLoQLM5Vuquijkm6UMFNTy1cMSiHLVoGIrIy4Ss8cHB30PX1z2GePxtbWVk488UROO+20zft6e3u5/PLLueWWW1iwYAF///vfeeCBBzZvh8NhQoXMEe1FMwH1wDjgt8DYhM8OuWi8JO2Mo79l8zXglwnb3wR+DjQAv8BYPTPT1DcDeBF4ccyYMWk1ejLy5T2SCbcv5xUVMQXzN5Vobs/B/fCdf85ofo96eSpX3asB6dPa6t6U5e65p01ho0KPwka955625JW5WMeypc3etG16wf1cfQnP7RUYN8NjbsnGsvGbqVOn6htvvOFbfUW3bFQ1CkSBM7JTbVmTLFiPqup64KJMX1bV2cBsMCkGvDYeCoWora3dbBWkexvwa5G2l8gc1dXilJOUSxRCoRBVVVPp6TmKQOD5lOcQDAa5/fZ/sHDheqalsTLcLOb2O7iBl3JOIIfM9fUG6FPo7k1dbv36CUAMqKCiIuBsJ6ksFsvYaDAIzU8FfF3I7y6liTvLdqgQXh6mu6+bPu2ju6+b8PKw79ZLIfHbUsoHXofRnlPVKSLSgXlL3nwI8/Df1lfpttAOjE7YHgWsylNbWxFfrBk3O9MNK02ZYp45tbW5xcT0MzKHUxKRZkAQUcw0WHK5rrhiAt3d8OyzMGFC9uGn3MrmVil5KVdb6299FRXdxGKVVFdXJi/ntlFg5syQ8zecsoxb3CvgMN3d3fT19WUOpeQCvyMVFTpCUmhciKpAFd193VQFqgiNC2VfmcUduZhF+fqw9TBaJbAUGM8WB4EDsqk7n95ofg4b+T1n7lY2L45h99zTpsce+1TyYSWP+L1soyjlXFbm5+Sw+yFD/7wW/b428+A86M5pMcPEv1tKYRjNb0rJG+0/gd1zaThN3b8FVmMiprcD5zv7vwS8hZmf+UG29ec3Npr/oaAKvQjey8PLzygIQ2k+we81Fm6vE79eDvz2VPf7RailRbW6plelok+ra/yZF0uHVTb5VTbXA68CzwKXADvnIkQhP9m7Pvv3Blys5Qx+noMbhwMvCikfK6pLkWKdq5/XXFEdOlyUu+i7yxXpMV7t0q0XfXe595P0QL6UTUdXh676eJWr8DJ+U3QHgYShtxuBG0VkInAa8LSItKvqMdnUV+p4SQXjNnOAr7lPXOLWq9VdZoPMThPGZdhMmqebJ/Dibl3uuO0T/9v1L5Oo36mW3Nbn+r4Z9zQEvgZ9CoEes830rYoVM5NGaG4IgPA54aTHO7s7eWv9W8Q0RkVnBfs07FM+qQRSkGvytLXAB8B6YKfcxSlN3Ho3eaHclzO4cZoIhUJUVFQQi8XSevEZT7kqV95+5Y7bPvG/XX8zifqZasltfW7LTT9xb+a0fomed49i2J7PM/3EWVuVce/FV5ycgx1dHcTUvJTENEZHV8fQVDbO4s7TMArm98C3VLUgCyyLgQdHoyFFpkWnwWCQ5557LqMXn1tvv2LidgFgJoLBIBMnTiQajTJ//vyCnatbq8DPt/iiWfCjg4SvmeX8XrOS/l5+u9v7zYjqEVR0VhjLRioYkS7lc7mQzdgbMAs4KJfxu2J98j1nYxl8+LkAUNX/UCV+UQ6h0fzyIPNzniirQJy/mKTjbx+f9jyymbO588479YADDtAxY8boHXfc4VmuOEWfs4mvswEuAy6Jj7OT/3U2Rafch70s2RNeHubTXifd9yBYAJiKYr3FuyWyMkLjvMbNa2Oapzdn/Tv4Pk/kgcjKCG1r2ohpjMZ5jSnPw2vK54ULF/LEE0+wePFiPvzwQyZMmMDFF19MZWWusyX+4DWCwBTnb3kPHlosHgiNC1FbWTvoFwD6Pa/jN36v+vdznsgL4eXhzfMxfr683HHHHdx7770MGzaMXXfdlWHDhm12RikFskoLLSI/dLPPYhkMBEcHaZ7ezM2fvzmnt+lSp9SzLsdX/QckUNZKPzQuRIWYR69f59HT00NbWxv77LMPAKtXr2bHHXdk/vz5PPbYY6gq5513HnfffXe/7U/jnk8FIFv76ovA9wbsOyHJPovFVyZPnlzwyXXIT+TfUqSUh4vjSt8PR41iEhwdZOLOE4luijL/lPm+nMdrr71GNBpl6dKljBs3jpkzZ3LZZZdx9NFHM2fOHN5//31OO+009tprr37btbW1PpyRO0ouxYDFkopIJEJbWxuxWIzGxsa0wSTjLsXlEKDQUtxw/8Wgvrqe+up6385l8eLFnHnmmZxxxhls3LiRU045hRkzZmw+Fo1GueCCC5JuFwqvls0DwGMYb7SrE/Z3qOq/fZPKYklCPhZEWqVUfPyc+B+qJMtnE6eyspLrrrsu5Xah8DRno6pRVV2uqmeo6nsJH6toLHknviASKOvFn9Fto6wYs4LIykixRSkJkk38D3bC54RTRg/IhtbWViZNmtRvXzQa5dJLL+Xss89mp5122mq70GS7qPPrwF9VtUNErsGkcb5FVV/2VTqLJYFiLYj0k8jKCG2T2ohJerfXwYKb4TEb7j93klnm9fX1/PznP0+5XWiydRC4VlV/LyJTgOOAHwF3A0f4JpnFkoTFixcXW4SUuHmwhpeHiUkMKjK7vfoVtaBYRFZGmPKrKcQ0Rm1lbUrFOlgm/i3pyVbZ9Dl/vwzcrap/FJEb/BHJYsmdaDRKNBolEokUxAJy+2ANjQtRoQmx0VK8xZf6PEYkEskYXii8PIw4SXYzKdbBMvFvSU1W62yA90VkNiY+2qMiUp1DXRaLr8S91pYtW0ZjYyORSP7nRpIt1EtGcHSQia0TGb9sfFoFEo9aUOh5jNmPzea4W45j9mOzU5aJp5i+9tpr0/bvYFkXY/GHbC2brwPHA/+tqhtEZBfgv/wTy2LJnmKE8Y8v1ItpLOODtf7jeuo/Tu/2WoyoBbMfm82Fz18IAXj8+ccBmHHCjK3KhcNhukZ2ERsTo2tFV8r+tcNjlkSyjY22FpM4K56HRJztQRsbzVI+FCOMv98L9YrxoF740kIIYMYo1GwnUzYNkxqIfTMGAYj1xWiY1JCyTjs8ZomTbWy0QRDv2jJYKZbX2uIL/XVeKPSDetoh04xFo0DMbCdjfd16KoZVECNGRUUF6+vWF0xGS/lSGuFALRafqa+vp76+viTdo0t1AWnciln40kKmHTItqVUDZoivurLauirnQtzaLtFrIR/4sc7mWmAydp2N75S766ul/JhxwoyUSiaOnYuxZIMf62yOZQissyn0g7/UXV8tpUMxXkrsXIzFK9m6K2+1zgao8kek0iO+huL7T36fxnmNBQkzMhRDeFi8U4xr0wuRSIRZs2YVxP28rIhGYcUKk4fbR+666y4OPPBAxo4dy89+9jNf686VXNbZ3MMQWWfjdg0FmJt/1rOzMt70mcrZNQq54Tb+WDQaZcWKFWX7MPRybRaaSCTClClT+P73v1+w9U5lQSQCbW2wbBk0NuakcDo7O1m9ejWdnZ39MnUuWrSIG2+8kd7eXh8Fz41sh9FOxayz+ZGzzmZX4Cr/xCot3K6hcLuK3E05Oy6ePZGVEVoPbgVIG38sEonQ+q9W6CNjyoJSpZSziIbD4fjSiIKtdyoLwmGIZ9DMIf92Z2cnb731FrGY8Qr8yU9+wpw5cwZXpk7gU2A4cIazPQzY4IdAXhCRPUTkPhF5MJ/tBEcHee7c52j6QlPGVd9u3jK9rDaf+bmZVtF4ZHN/SobfIRw2A8K65WFYbpRyFtFQKERVVRWBQKCso3T7TigETvTyXPJvd3R09Fu8/K9//WurTJ3RaJRzzz2X9vZ2zjvvPO69996yy9R5FxADvgDcBHQAC4HD3FYgInOAE4G1qnpgwv7jgZ9ilpf9UlVvS1WHqi4Fzs+3sgF3E6JuLSC/30at11p/NvenQlVlmt8hFDJXWV95pyxwO1kfmhsC8DW0fTqCwSDNzc0ZY6gNOYJBmDjRzNvMn591atQRI0ZsXry8fPlyPv74460ydY4cOZIxY8Zw5ZVXct9997FmzZryyNSZwBGqerCILAZQ1Y9ExKuDwFzg58C8+A4RCQB3YtJOtwMviMgjmEfCrAHfP09V12Ypf16IW0CZHvx+DpG5HbobLLhRrMHRQSa9PInodlHm35p6NX8wGGT4KcPpbe+l+UflN4TmlWhXlOimKJGVkcJ5rQWDg75fs6K+3nxy6Ju6ujr22WcfOjo6ePHFF5Nm6uzs7GTp0qVUVlZSV1dHXV1d2WTqjNPjKIZ4yJqRGEvHNar6jIiMG7D7cOAdx2JBRH4HnKSqszBWUMnj9i3TL9fRZENyQz1kPbiLPwZQuWsllbtWDvoHYmRlhLY1bcR0aOTRGSrEFcgbb7yxVabO3t5eLr/8cm655RYWLFiw2cIsi0ydCdwB/AHYSURuBZ4DmnyQZ3dgZcJ2u7MvKSLSICK/ACaLyMw05WaIyIsi8uK6det8ELN0iA/JZfJaK3UXWTeUsudVqWP7rsQIh32NHpAsU2dlZSVz5sxh7NixXHzxxey6665ceOGF5ZOpU4x7yTPAS0AjJgjnyar6ug/ySJJ9mqqwqq4HLspUqarOBmYDHHrooSnrK0fcDsl5sYD8xM95gnx4XvUGeukd1lvQoaVi4CUqdTFwkx/Hkpp0zi2JXmsXXnjhZieCQuNZ2aiqisjDqnoI8IbP8rQDoxO2RwGrfG5j0OHWeaEYLrJ+zhP47Q4eWRlh44iNQHoX6cGA31Gp/SSeH6e7u5uqqqqydEHPB52dnXR0dDBixAjq6uqyrifRay0Wi9HR0ZFTfdmS7ZzNIhE5TFVf8FUaeAHYW0TGA+8DpwPf8LmNIUk+1u1kmgPKxzyBn2FSkrlIl9JD2G/qq+upr848j1VowuEw3d3d9PX12fU4DgPX0Oyzzz5ZK4hEr7WKigpGjChO0P5slc3ngQtF5D1gI04+G1Wd6LYCEfktEAJ2FJF24HpVvU9ELgX+hvFAm6Oqr2Ypo2UAbh/UbhwJ3EzWF2vozi1uXaQHC4VyefZKfD1O3LIpVxd0P/HTGkn0WsvVSsqFbJXNCbk2rKpnpNj/KPBorvVbssOtx5cbRVLq8wTB0UGGdwynd1gvzecP3iG0Useux9kav62RuNdaMclK2ajqe34LYikNwsvDiOOnkc4acTMHVMrzBHEq+yqp7KssSdmGEnY9Tn9KxRrxE5s8zdKPeADQTI4EbueASnWewFKeDCWvtVKwRvzEKhtLP7w4Evg5WT/5nsklbQFZik+pe6259R7r7Ozk2GOHEQhU8uyzgQJKWFyyzdT5Q1X9XqZ9lvLEV48vF5PSdnW7xQ3hcHhz4MhS81pz6z0WL9fVtTci3XR2yqCyXtKRbQSBLybZl7PTgGVoYle3W9wQCoWora0tySjSybzHMpVTJWW5wYgnZSMiF4vIEmBfEWlL+CwDluRHRMtgJ+61BpSk15qlNIh7rd18880lN4QW9x4D0nqPxct1dlbwwQdVvPZava9ylHKmTq/DaA8Aj2EiMF+dsL9DVf/tm1SWIYXbaNkWS6l6rbn1Hqurq+Ojjz7DO+/UEovBV74iNDfnFPx5M4mZOj/88EMmTJjAxRdfTGVlaUzNe5JCVaNAVETOBU4BxsXrEBFU9SbfJbQMCfycJ7JYioFb77F//nMbPxJ1bsUdd9zBvffeO+gydT4MnAT0YiIIxD8Wi8UyqIhEIsyaNYtIxJ9I6T4l6uxHT08PbW1tW2XqnD9/fr/MnCtXruyXubOnpyf3xl2SrX01SlWP91USi8ViKTHy4W7tU6LOfrz22mtEo9GtMnUeffTR/TJzjh49ul/mzmHDhuXeuEuytWxaRGSCr5JYLBZLiZEsSKgf1NfDmDH+KBqAxYsXb87UOXHiRMaMGcOMGTPYc889Wbx4Ma2trRx33HFbZe4sJNlaNlOAc0VkKdBFFoE4LRaLJV/4FWmgXIKEtra2bpWpM048M2eqzJ2FomiBOC0WiyUfRCIRpkyZQiwWo7a2Nqehr3wFCfUxSSdglM3FF1/cb180GuUHP/hBv8ycc+bMAeCqq67yVwAXZKtsVgBnAnuo6k0iMgbYBbABOi0WS1EJh8OYhML+RBooVXfrRJIN79XX1/Pzn/+88MKkINs5m7uAIBBPE9AB3OmLRBaLxZID8aGvUow0MJTJ1rI5QlUPFpHFAKr6kYhU+SiXxWKxZIXNj1OaZKtsekQkACiAiIwESmf1kMViGdKUw9DXUCPbYbQ7gD8AO4nIrcBzQJNvUlksFksZoarFFsE38nUu2WbqnC8iLwGNGLfnk1X1dV8ls1gsljzil3t0TU0N69evp6GhYbNjQrmiqqxfv56amhrf6/asbMT05ihVfQN4w3eJLBaLJc/4GRlg1KhRtLe3s27dOp+lLA41NTWMGjXK93o9KxtVVRF5GDjEd2ksFoulACSLDJCtshk2bBjjx4/3WcLBR7ZzNotE5DBfJbFYLJYCYd2jC0+23mifBy4Ukfcw0Z5tuBqLxVI2WPfowmPD1VgsliGJdY8uLJKtm5uIbA/sDWx2W1DVZ3ySK2+ISBR4O2FXPRB1ub0j8KGP4gxsy4/yqcok2+9mX+J2PvsilTy5lE93PJv+GErXRrL9g+le8fvaGLhdTn2Rrkzi/rGqOtJDu/1RVc8f4AJgCfAR8BTwKfBkNnUV+gPMznYbeDGfsvhRPlWZZPvd7Btw/nnri3z0R7rj2fTHULo2Mp1vufeH39dGumul1Psim2sjm0+2DgLfAQ4D3lPVzwOTgXLx+/tTjtv5lMWP8qnKJNvvZt+f0hzzG7/7I93xbPpjKF0byfYPpnvF72tj4HY59UW6Mr6dR1bDaCLygqoeJiKtmDhpXSLSqqqT/BKsFBGRF1X10GLLUQrYvuiP7Y/+2P7Ygu0LQ7YOAu0ish3wMPCEiHwErPJLqBJmdrEFKCFsX/TH9kd/bH9swfYFOTgIbK5AZCpmEumvqtrti1QWi8ViGVRkZdmISA3wbUx6aMUE4sx2/sdisVgsg5xs52wWYBKm/cbZdQawvap+3UfZLBaLxTJIyFbZvKKqB2XaZ7FYLBYLZD/0tVhEjoxviMgRwPP+iGSxWCyWwUa2ls3rwL7ACmfXGOB1TLZO1SEYI01E9sOsP9oRaFbVu4ssUtEQkZOBLwM7AXeq6uPFlai4iMgewA+AelX9WrHlKTQiMhy4C+gGwqo6v8giFZWhej1ka9kcD4wHpjqf8cCXgBOBr/gjWuEQkTkislZE/jVg//Ei8qaIvCMiV6erQ1VfV9WLgFOBsvWp96kvHlbVbwHnAKflUdy841N/LFXV8/MraWHx2C+nAA8618R/FFzYAuClPwbj9eCGbJXNGmAa8BPgx5iLaY2qvqeq7/klXAGZi1GgmxGRAHAnJujo/sAZIrK/iEwQkT8P+OzkfOc/MJ55zYUV31fm4kNfOFzjfK+cmYt//TGYmIvLfgFGASudYn0FlLGQzMV9fwxJsl3UOQ/jjfYzZ/sM4NdAWXqjqeozIjJuwO7DgXdUdSmAiPwOOElVZ2EsuGT1PAI8IiJ/AR7Io8h5w4++cLK53gY8pqov51nkvOLXtTHY8NIvQDtG4bQySJdIeOyP1wosXkmQ7Q+/r6qer6pPOZ8ZwD5+ClYC7M6WtzEwN8zuqQqLSEhE7hCRe4BH8y1cgfHUF8BlwDHA10TkonwKViS8XhsNIvILYLKIzMy3cEUkVb88BEwTkbvJf3y9UiJpfwyh66Ef2Vo2i0XkSFVdBIPWG02S7EvpTaGqYSCcL2GKjNe+uAO4I3/iFB2v/bEeGIxKdyBJ+0VVNwLnFlqYEiBVfwyV66Ef2SqbI4DpItLPG01EljB4vNHagdEJ26MYGvHfkmH7oj+2P5Jj+6U/tj8SyFbZHJ+5SNnzArC3iIwH3gdOB75RXJGKhu2L/tj+SI7tl/7Y/kggqzkbx+PsY2BnYGz8U67eaCLyWyAC7Csi7SJyvqr2ApcCf8OsIVqgqq8WU85CYPuiP7Y/kmP7pT+2PzKT7aLOCzALGOMeJkcCEVX9gq/SWSwWi2VQkO0wWjxT5yJV/byIfAa40T+x8seOO+6o48aNK7YYvvFB5we8//H7ZkNg9xG7s0vdLgVp+8033wRg3333TVmmdWUrAJNGTypYm17KWSwWd7z00ksfqurIbL+frbLZpKqbRAQRqVbVN0SkLO7qcePG8eKLLxZbDN+IrIww5VdTiGmM2spafj/99wRHBwvSdigUAiAcDqcss90V2wHw4u3+9LmbNr2Us1gs7hCRnKZIbKbOMic4OsjEnScS3RRl/inzC6Zo/MYqB4tlcJOtg8BXVXWDqt4AXAvcB5zso1wWD9RX1zOmfkzBFU00GmXFihVEIpFB3SZAaG6I0NxQQdu0WAYTOYeOUNWnVfURmxJ6aBGJRGhra2PZsmU0NjamfPj3ru5l0wubfFEObtu0WCylR7bDaIOKnp4e2tvb2bRpU7FFyYrrD7gegNdffx2AmpoaRo0axbBhw/qVi7+Zh88J59xmOBwmFosB0N3dTTgcJhjsb1lFIhE2PrQR+qCxsZHm5uatysSJRqNEo1EikUjKMm7atFgspYlnZePkptikqoMmemt7ezsjRoxg3LhxmBiS5cV+7Lf5f1Vl/fr1tLe3M378+Ly1GQqFqKioIBaLUVVVtXnOJZFwOGxi/Gp65RC3WGKxWFql5KbNOK2TWrM+t2REu6JEN0WJrIyU7byYxVJMMg6jiUiFiHxDRP4iImuBN4DVIvKqiPyPiOydfzHzy6ZNm2hoaChLRTMQEaGhoSHvVlowGGTixImMHz8+rXJgDPA5CIwLpFQOySyWbNvMB5GVEdrWtLFswzIa5zUSWWmH7ywWr7ixbJ4C/g7MBP6lqjEAEdkB+Dxwm4j8QVV/kz8x889gUDRxCnUu9fX11NfXp37ojwLOdmSqFLOdBC8Wy7KpywAKOnwWXh4mpo4y7OsmvDxsrRuLxSNuHASOUdWbVbUtrmgAVPXfqrpQVacB/5c/ES1+Ee2KsiK6omBv5uHlYRP3tgJ6tddsJ6FYFotbQuNCVIi5VaoCVYTGhYorkMVShmRUNqraAyAiXxeREc7/14rIQyJycGIZS+lSjKGgzQ9lzfyQrq+vZ8yYMSWnaGDLWqbx242neXqztWoslizw4vp8rap2iMgU4FjgfuDu/IhV+oRCobTDPaVGsqGgVLhdU9I6qTXtRHxwdJDhHcOp3lTt20O6N9DLpppNGZWl23JuKdZaJotlsOBF2cS9z74M3K2qfwSq/BfJ8tBDD3H55Zf7Wmc+hoImTZrEpEmT0pap7KukZlNNxod0OBzOGD0gsjLCxhEb6arpSmuduS1nsVgKhxdl876T8vg04FERqfb4fYtLFi9ezMEHH+xrnV6Ggvyc23GjkNyy2RqT9NaZ23Ke2j4n7Mv6JCg/q9hi8QMv62xOxSRN+5GqbhCRXYGr8iNW6eNmEaJX3nrrLS655BIWLVpEQ0MDGzZs4IorrvClbjBDQfXV9WkVTXxuJ6YxGuc1ltQcRb85oMrU1pnbchaLpXB4sUw+BYYDZzjbw4ANfgtUDuQjbEpXVxennnoqP/7xjxk5ciSLFi3ipptuKnhUAy9zO4X2bguODjJp10mM3z69dZaPuSKLxZIbXpTNXZgkaXFl0wHc6btEZYDbRYheeOKJJzjooIPYbbfd2Hbbbdlll12oqamhr8+/QA1uhoLczu0Ua6Gj24l6t3NFFoulMHhRNkeo6iXAJgBV/Ygh6iAQX4QIZFyE6JbW1lYmTJjAK6+8wsSJE1m7di0jRozgk08+4dxzz6W9vZ3zzjuPdevW9dvu6fHX69zt3I4XC8hisVi8zNn0iEgAUAARGQnE0n9lcBJfhBiNRpk/f74vczbbbrstbW1tVFZWMnHiRK6//nouueQSRo4cyZgxY7jyyiu57777qKur67c9MNimH7iZ24lbQDGN+eLd5jZIqF+T9BaLpbB4UTZ3AH8AdhKRW4GvAdf4KYyIzAFOBNaq6oHOvh0wEQrGAcuBUx2rChGZCZyPccu+XFX/5qc86cgYqsUjZ511Fl/96ld56KGH2H777Tn99NO57LLL6OzsZOnSpVRWVlJXV7fVdj5w80AfLEnbLBZLYXA9jKaq84HvArOA1cDJqvp7n+WZi/F4S+RqoFlV9waanW1EZH/gdOAA5zt3OZZXWbLDDjvw9NNPM3r0aP75z39y66230tfXx+WXX84tt9zCpEmT+Pvf/95vu9hZLe1CR4vF4hZPKQZU9Q1M1Oe8oKrPiMi4AbtPAkLO//cDYeB7zv7fqWoXsExE3gEOBwoyU52PB31XVxcdHR00NDQAUFlZyZw5cwC46irjZX7MMcf02y51ijXsNal1UlHatVgsyXFt2YjhLBG5ztkeIyKH50+0zeysqqsBnL87Oft3B1YmlGt39m2FiMwQkRdF5MV169blVdhcqK6uZunSpcUWY0hhF1haLIXBi2VzF8Yh4AvATRjX54XAYXmQyw3J4uhrsoKqOhuYDXDooYcmLWPxjp2st1gsbikH1+c1TrQCnL9rnf3twOiEcqOAVQWQx1IGRKNRVqxY4cuCW4vFkjtelE2xXJ8fYXMKLs4G/piw/3QRqRaR8cDewD8LII8lD/gZjSAfER4sFktueFE2A12fnwOa/BRGRH6LmeDfV0TaReR84DbgiyLyNvBFZxtVfRVYALwG/BW4RFX9W25vKRh+RyPIR4QHi8WSG67mbMTkGX4GeAloxMyXnKyqr/spjKqekeJQY4rytwK3+imDpfD4nXbZS5ppi8VSGFwpG1VVEXlYVQ8hj67Ppcybb5q/++5bXDkGI35HI8hHhIehQiQC4TCEQmC7zeInXrzRFonIYar6Qt6kKWH6+synsxPq6tyHV0nFa6+Z+saPN/UNZfIRjcDvCA9DgUgEpkyBWAxqa6G52Soci394mbP5PEbhvCsibSKyRETa8iVYKdHZCZ98Al1d8NZbZjuf9eUjU2epM5SiEZSqp1w4bBQNQHe32bZY/MKLZXNC3qQocTo6tvwfi/XfHogbiyVZfYll85Gp01IaxD3lYrEYjY2NNDc3l4z1FQoZi6a7G6qqzLbF4hdelM20JPuiIvKSqrb6JE9JMmLElv8rKsx2tCtKdFOUyMrI5rfxuMUCxmLZZ5/kCidZfeY7+c3UaSk+yTzlSkXZBINm6MzO2VjygZdhtEOBizAhYXYHZmBilt0rIt/1X7TSoa4OttkGqquNAlnyUXJXXbcW0MD66upSZ+pctmxTxmG7zk5YvTr34T1L/slHLiQ/CQZh5kyraEqB2Q8v4bgLw8x+eEmxRfEFL5ZNA3CwqnYCiMj1wIPA0RiX6P/2X7zSIRAwn7o6CC9O7qqbymJJxv77998emKmzrm4XAoEa2tv7+Oij1FZSZye88caWNlOVywel7LnU2nq7r+X8Ih+ecpFIhHA4TCgUKhkryZIbsx9ewoWn7gm9+/H4r7phwRJmnDyh2GLlhBfLZgzQnbDdA4xV1U+BLl+lKkH23XeL23Oq1Ml1dfCZz8Duu3t/6A/M1Lls2Vq22WYETzyxgGeffYyPP1bOO+887r77bh577DFUzfbatZ9uriPTfNLkybDHHkZJpMNNubjn0ve/D42NmevMRPSnYVbcFM5YTyQCs2Zlbq+3dxs2bdrJt3J+Ul9fz5gxY3xTNFOmTOH73/9+SUZLcPt7DQYikQizZs3y5TdY+Nh66K0CrYTeYWa7zPFi2TyA8UaLh4v5CvBbERmOWcU/ZEjnqltXl51lMTBT509/ej2nnnoJBx98NH/60xwCgfc57bTT2GuvvZgzZw7vv2+2d9qpln//29SRzpqKRKCtzSikxsbUbq1uyyXzXEpWbvJkiEZh/vzU1o8X2dy45kYisHHjXoBkrM9NuVKmlOeAIhHze8YdDvzoX7fWdKGt7rjSj8Vi1NbW5uz4Me2EBmPR9CpU9jDthAYfpS0OXpKn3Qx8C9gARIGLVPUmVd2oqmfmSb6SxW9X3bPOOou3336bm266ibvvvpudd96B6667jCOO2JP29sW88UYrxx13HHvuuSeLFy+mtdVsu7Wm3Lq1ui0XChnlBqk9l+JKZNmy9NaP37KZ/QKIL+VKmVAoRG1tLYFAoKBzQG4slnAYPv3UeGf60b9x5XXttemvJ7fl/CQcDmMCrfgTImnGyRO4Z8G7HPut57lnwbsph9DcWo6lYGG6tmyckDX7AfWqelM8n42q2uCXPhDP1HnQQQfx5JNPbk6gNmIEVFdXct11120uW1nZf9uNNRVXDrFYerdWt+6vwSA891z6t0e31o/bNr2cQzzbRFWV5FyuWLiZiwkGgzQ3Nxd0zsatxeLFldqNJRIOm7oSlVeq685NOT8JhUJUVVXR3d2dUelHVkYILw8TGhdK+7I64+QJzDg5dZuRCHz+C310dwtVVcpTTwZSWvB+W5hZoaquPsDdwJ3A68729sALbr9fKp9DDjlEB/Laa69tta8YbNq0ScePH795e8OGDXrJJZfogw8+mHQ7HcnOqaVFtanJ/E2H23KZaGlRrahQBdXa2vT1+S3bpEkX6/jx92QsN3z4m1pd/X7GclOnmo8fTJ06VadmqKylpUUrKioU0NraWm1JI6Bfv5eqasuKFm16pklbVqSurKnJ/KagGgiY7Vxka2kx10cgkP46yUc5v/rN1NeiTU1N6X+rFS1ae0utBm4MaO0ttWn7ORMXfXe5Ij3mt5Buvei7y5OWa2oyfeHm90oH8KLm8Oz1omxedv4uTtj3Si6NF+NTysrGT0rlnPy+od3i5oHurZx/yqa+frHW1y9OW6apqUkxZpcGAgFtSvGEcPtgdUPLihatuLFCuYG0D0I/21T19jD066XE73NwS9MzTRq4MaDcgAZuDGjTM1k++VX1orvuVyo3KtKtVG7Ui+66P2m5lhbV6ppelYpera7pzfpcc1U2XhwEipXPxlLGBIPlN+leCphhman09BxFIPB8ymGZ+LwIZB4yyjRU5Tb6tt+LP0MhM7zjZrjN7fWUqVwxhtrAeLJWBaro7uvOOejs9BP3Zk7rl+h59yiG7fk800+clbzgqAg6fSa8exS65/MwahZQ+JvSi7KJ57PZ2cln8zXgmrxIZbEMeYL09DSjKogoEEhayu28iBtPPi/Rt/18iShG5AIvCs7PdUzB0UFuP+AfLHxsPdNOaCA4Ovu1M8HRQcLXzHLmf2alnP8JLw/Tt/tz6G5P0yeBnFN4ZItrZaOq80Ukns8G8pDPxmIpVaJR84lE0lkO7h5Kvb3b0Ntbl7aucBhUA0751G/ebh/Ubpw1gqODPHfuc64mr/2m0Baw237z26U5EoErvjGB7m549tcwIcfJ+uDoYMbfyU9rKhcyKhsR+X8pDp0gIieo6o99lsliKSncrAOKRCJ89rNXA0dRVTWTcHhW0oeS27U9Xjy53Dyo3XryuXl4DRbc9Jvf65iKMXwXHB2keXpzUV4iEnFj2cSXCe4LHAY84mx/BZO902IZ1LixCubNext4DDDur/PmPZj0oZRsbU8uFotb3LiqW7Ymvo7JjUuzu/rcD9/5SSm8RGRUNqp6I4CIPI6JjdbhbN8A/D6v0rlARI4HfooZ1P6lqt5WZJEsgwx3VsFUoApzS6mznbwut2t7/B5amvlmCHaFmcGwf5UOcvxexzSUI2t7cRAYGButGxjnqzQecbzj7gS+CLQDL4jII6o6pMLnWLKndVJrxjLBIEycmD7szvTpY/nVr/ro7u6jqqqC6dPHpqxr+PB36O2to7l5tyH1sClXgsGgr4tlh6qHphdl82vgnyLyB8yr2VeB+/MilXsOB95R1aUAIvI74CSGWKw2S/6przefVA+JYBCeeirg6o21svITKis/IRjcLR+iWiwliRdvtFtF5DHgc86uc1V1cX7Ecs3uwMqE7XbgiIK0HB//KMeAWpa8MFTfWC0WN7jxRhNn9Siq+jLwcroyBUaS7NtKDhGZgUn2xpgxY/ItU8489NBDhMNh7rjjjmKLYnFobW11/puUc12TJl3h/BfOuS6LpVxwY9k8JSILgT+q6or4ThGpAqYAZwNPAXPzImF62oHRCdujgFUDC6nqbGA2wKGHHloMpeiJxYsXc/DBBxdbDEsCg0FBbFGYFkvhcaNsjgfOw+SuGY9JMVCLSU/wOPATVW3Nl4AZeAHY25HrfeB04BsFadnNKj+PvPXWW1xyySUsWrSIhoYGNmzYwBVXXOFL3ZbSIdfw89nSu7qX3vZeIpFIyeS8sQwdMuazUdVNqnqXqh4FjMVEEJisqmNV9VtFVDSoai9wKfA34HVggaq+mveG3SZq8UBXVxennnoqP/7xjxk5ciSLFi3ipptuYtOmTT4IbBnqRCIRNj60ka5IV0lm9BwMRFZGmPXsLCIr/elbPzN/lgJevNFQ1R5gdZ5kyQpVfRR4tKCNuk3U4oEnnniCgw46iN12241tt92WXXbZhZqaGvr6+nIW12IJh8PQB2jpZfQcDERWRmic17g5JEzz9OacFlH6HSanFHCdqdOSgJs0lR5pbW1lwoQJvPLKK0ycOJG1a9cyYsQIPvnkE84991za29s577zzuPfee3nsscdQVc477zw+jYf8tVjSEAqFzLJnoaAZPYcK4eVhuvu66dO+zRGzc6ovSZiccseTZWNxcLPKzyPbbrstbW1tVFZWMnHiRK6//nouueQSRo4cyZgxY7jyyiu57777WLNmDXPmzOH999/ntNNOo7a21ocTsgx2gsEgNWfW0NPbw+1n3l72b8l+4Gc0Z7+DXfodJqcU8KRsHA+0gKra1+lMq/w8ctZZZ/HVr36Vhx56iO23357TTz+dyy67jM7OTpYuXUplZSV1dXXU1dWxePFiotEoF1xwgS9tWwY/kZURuvfoJqYxrnj5CiYcOCHnYZ5CpqL2G7+HqfwOdlmMdN/5xrWyEZHvANcBm0TkY+BOVf153iQbYuywww48/fTTHHTQQTz55JM0NDTQ29vL5Zdfzi233MKCBQs2X3iVlZVcd911xRbZUka4TYzmhsEwn+B3NGfwP9il32Fyik3GORsRuV1EpgPfAfZT1d2Bo4H9ReTmfAtYsoTDvkcP6OrqoqOjg4aGBgAqKyuZM2cOY8eO5aqrrmLy5MlceumlnH322ey0006+tm1JTzQaZcWKFWXrGRRPjAbkPMzjZT7Bbw8tvwiFQlQ4866DZZiq1HFj2TwNTAZ2BFocq6YNWAJcJCL/q6ob8ifi0KG6upqlS5emPF5fX8/Pf26NyUITiURoa2sjFovR2NhYlm/ywdFBJu48keimKPNPmZ/TG7jb+QQvHlqRlZGC5lsJBoM899xzroapMqXTtrjDTYqBPwB/EJEjgf/EuD4fBEwEdgDCIlKnqnvlVVKLpUjkY8ilGNRX11NfXV+w+YTw8jCf9prp3XRDd5GVEab8agoxjVFbWZuz27Bb3AxTRSJmKV08/0yqZHeWzHhxELgEWAC0Yqya/YAlqhpyHAcslkFJfMglFosN+iGX+LllcrV186AOjQtRW1mb0UMrvDyMOGEOc51P8ptwGOKrCwqVWXOw4nqdjaq+jYmo/CAmXE0bJs0Aqtqd5qsWS1kTDAaZOHEi48ePL8shtGIRHB3k9oNvp7GikdsPvj2lAom7DQckkHE+qdBzQPH03IFA5iV1kQjMmuVLQJFBidcIAt3AX5yPxTJkqK+vp76+ftArmmg0SjQa9SV+WiQS4eITLyYWi/Fs7bNMaJ6QtE63bsN+r9J3g9vMmna4LTM2goDFYgG2OEIsW7YsY/w0NxaGF6+14OggMz83M63yiM8BuVml76cFFAzCzJnplUc4bBRNX9+W4bZCUE7x02wEAYvFArh3hHA7oe/3XJfbOaDIygifn/t5umPdVFVU8dQ5T6V0THDjAeemXChkLJq4ZVOIab1IJEJjY+Nmr8BSH+L1sqhTgDOBPVT1JhEZA+yiqv/Mm3QWi6VguFUObheIenEvdkN8DmjhSwuZdvC0lA/+ec/Mo6u3Cyqgq6eLec/MI3hm/7Juh+TclnM73OYn4XCY7u5u+vr6ysJL0otlcxcQA74A3AR0AAuBw/IgV8ljs0JbBhtxR4hoNMr8+fNTPrjiC0RjGss4oe/nKni3c0AsZ3OEa2LO9gCSBc5MpkTcloPCpwUPhUJUVVWVTfw0L8rmCFU9WEQWA6jqR9blOT/YtNCFYyglFAufE85Yxo0jRHB0kOfOfa6gizDB/TDf9C9MZ84359Czew/D3h/G9F9P36qM28CZfgfY9JNyi5/mRdn0iEgA876AiIzEvDdYfMamhS4M8YRi9OFPZIBimbuTJ/sWgdxtKHu/44BNnjw5s0XlMnJBMBgk/Otw2oewWw84vwNs+k05xU/zomzuAP4A7CQitwJfA67Ji1RlQB6yQtu00AXGS0KxouQTcRMnJZ41NhYzvrdl6HPrNhyQlzd5Nw9htwrTb8U6VHGtbFR1voi8hEkLLcDJqvp63iQrYfJxf8fTQv/617/mpJNOoqWlhf3335+LLrqImpoafwS39GNzQrG+EgzGGInAlCnmIqutTX2R5SFrbKHxEg6onN7kLf3xtM5GVd9Q1TtV9edDVdFA8vs7V2xa6MITDAYZfspwqoPV/riNRqOwYoU/S8jdXmRelriXKPHhsUAgUHpK3+IbrpWNiHxdREY4/18jIg+JiC8TC07dr4pITEQOHXBspoi8IyJvishxCfsPEZElzrE7HNfsgpCHrNAp00IvWLCgXxrolStX9ksT3dPTk3vjQ5jKXSupOawmd0UTN3eXLTPmbq4Kx+1FFve5vfnmshxCgy3DYzfffHPJrxWxZI+XOZtrVfX3IjIFOA74EXA3Jl5arvwLOAW4J3GniOwPnA4cAOwG/F1E9lHVPqftGcAi4FHgeOAxH2TJSB6yQqdMC3300Uf3SwM9evTofmmihw0blnvjg5CCz7H4PZwVDMJzz7lbuFFon9s8YIfHBj9elE18POfLwN2q+kcRucEPIeJDckmMk5OA36lqF7BMRN4BDheR5cC2qhpxvjcPOJkCKRvwPSt0yrTQItIvDfTANNGWEiFuicRi/pm7g0CJWEqDUkjj7UXZvC8i9wDHAD8UkWryH1ttd4zlEqfd2dfj/D9wf1JEZAbGCmLMmDH+S+kDydJCx4mngU6VJtpSAuTD3LVYfKBU0nh7UTanYoaqfqSqG0RkV+Aqt18Wkb8DuyQ59ANV/WOqryXZp2n2J0VVZwOzAQ499NCU5byQj1GagWmho9EoP/jBD/qlgZ4zZw4AV13luusthcJvc9cyKCm0lREOhzePGhUzrI0X1+dPROQpYG8ROdrZvcnD94/xKhzGYhmdsD0KWOXsH5Vkf1kzMC20TQNdXrR+0ArApKJKUaYMkdzLxbAySiWsjZdAnBcA38E82FuBI4EIJlZavngEeEBEfoxxENgb+Keq9olIh5Oq+h/AdOBneZTDYrFkg9uFqW6TwbhVSiWqvIqRYrxUwtp4GUb7Dibo5iJV/byIfAa40Q8hROSrGGUxEviLiLSq6nGq+qqILABeA3qBSxxPNICLgbmYrKGPUUDnAIslGVdcPQmAcFGlSIOPYW1c4VaJuM297LY+L+X8VFwuyrkNueM3peDt50XZbFLVTSKCiFSr6hsisq8fQqjqHzChcJIduxW4Ncn+F4ED/WjfYilnQnNDQIZAm27DXvj5AE6WUSxZ2fjC1EzJYNzW56ZcPhSXi3LBYJB/3H476xcupGHaNCYUyjorAUvPi7JpF5HtgIeBJ0TkIwbBPEkcVU3mer2ZNz98E4B9d/RFv+YVVV98ICxFxpUScYubdUBuQ+S4fQC7zSjmNhmM2/rclPNTcXkpF4kw4YorTJlnn4UJE/JvnZVIzmovDgJfdf69wXEUqGeQDF3V1NSwfv16GhoaUiqcTz75pMBSZYeqsn79ehtPzdIfN+uAwmGIX/+ZHqxuhr28ZBRzs6bIbX1uyvmpuLyU89s6cxs/z02becaLg0A1MA0Yl/C9SZhEamXNqFGjaG9vZ926dSnLrP33WgBe31j6IeFqamoYNWpU5oIWX4l2RYluihJZGUkfJdjl3Inb+lyVcxORwMuD1c2wV7xdPx9sbuvLVM5PxeWlnN/WmZuXg2LkrE6CuB1yEZG/AlHgJbZEE0BV/zc/ouWHQw89VF988UXP39vuiu0A2HD7Bn8FshQVv37XyMoIU341hZjGqK2sTZk+2O3bqNv6XLfr+kTK29urLPCrj/PhxZcGEXlJVQ/NXDI5XuZsRqnq8dk2VO50beiip7eH2Y/NZsYJM4otjqXECC8PE1PHpTVd+mCXMdTc1ue6Xbf4ZTlYUlNo68xLm3nES7iZFhGZkDdJSpjZj81m06hN9O3Zx4XPX8jsx2anLDt58mT22GMPIhmi/oZCIRtqZhARGheiQsztlDZ9sMtozm7rc92uZXASDMLMmUVXJG7IqGycMP5twBTgZSfUf1vC/kHPwpcWmiRbFeaz8KWFSctFIhFa/9XKsuXLaGxsTKtwxq5axZfb2lgyO3fF5Rar4PJHcHSQiTtPZPx249MPZcXnTpqa0g57uK3PdbsWS5FxM4x2Yt6lKHGmHTKNx59/3ERfi5ntZITDYY7shRDwbFdXytXBS2bP5u6336YK6L7wQpYAE2b0H5qLRCJs88orHK/KzFCIWelWGhdr/HyQrQPIlfrqeuqr6zM+8ENvzoRdIRycmbbc4gsX+9quxVJMMiobVX0PQERqgG9jLBwFnsPklBn0zDhhBg/87DKOWtXDYdOv4uQUczYnNjTwHTBKJBbj3YTIzYmsX7iQ/TCdr842A5TN2/Pm8YSqqau7mwfnzUuubNy6P2KsqQM//JAls2dvpdy2qtMv3303nlcezsHvCew/P9hp/rk9dVV+49przWIZRHhxEJgHdLAlBtkZwK+Br/stVMkRidD8124qFOSan0Hw5KQPsAnr1xPDjLYFKiqYsH590uoapk2j+/HHUUyuhIZpW1tKUzFKK66QpqaSLRxGYzETBjvNhLMbayp+rq599zOttXC7at1t4jEviw5dKq8RXUp9t5rvFMCiiqyM0LamjZjGaJzXaIe+LEMGLw4C+6rq+ar6lPOZAeyTL8FKinB4S16DDPngu6sq6K0Aqa5OOfk7YcYMTvxCDTeNCfDuPfckfeiPnT4dqa6mT4SK6mrGTp+ess0uEXqBvsrKlG2uX7hws/Ia5mynOtde58GvXV1pz5XaWggE0q8DGKhEsq3LS31uy0UiHLg+xtiP1Z9Uzi5I5j2Wc53nhP2JMmCx5BEvls1iETlSVRcBiMgRwPP5EavECIXoHVYBvTEqM4TduPCKPTnw1Q85/oLb0sY9in3zCJ4Fbj4nxXBWMEjgqacyDgVFgP8CjgYiqswCkpV0Y00BLGloYE+MNdXjDAUmdUF0u0rbzeI/Lwvi3NTnNmtmOGysVci8sjpeRyrF5ZK491hMY9Z7zDKk8KJsjgCmi8gKZ3sM8LqILAFUVSf6Ll2pEAxSE86cDz6yMsJvhr9L7LAY1796Bc2HTMhpiMTNRHI4HKZFlRYg0NeX0ilhwowZnP2jH3Hghx9y/G23pZyz+fP69TyC4+RQUcGX169Prmyg8OsAvKzmzrRaHsxLhECFAsMqCaTz1ItGzSfH4ba491h0U5T5p8y3Q2iWIYMXZTNkF3QCrh6Gfi+wczORHAqFqKioIBaLZQxZ/t5uu/HebrtxVRrngFAoxDUVFSyKxaitruZ/cnWVLtVwJUBkFFx9Nhy1Alr2UmaNSm4Vup17cjuUZb3HLEMRL4E438unIIMBL0MkmRSJ24nkYDDIxIkTiUajzJ8/P7V7NMYKyoSX+sqd8PIwz4yFZ8ZBQPpyXvXvul07v2IZgnixbCwZcDtE4kaReLGSFi92tx7DLX7XV6psfhlQqKp0seo/0xyQxWJJietAnIMFEYkCbyfsqscEGHWzvSPwYc5CjGAX6tjd2VI6WUUHH/QrU8VwGtgHM3+trOctuqkcIFsyBsqfbr+bfYnb/vdFZnlyKZ/uuDlWxXCqGUEXHSn6tx6IjoDhI2BkB6zrgI1J6ve7P/zui3Rl3O4v/L2SWcZsy2e+NrzvK9S9UqxrY6yqjvTQbn9UdUh9gNnZbgMv5lMWP8qnKpNsv5t9A84/b32Rj/5Idzyb/hhK10am8y33/vD72kh3rZR6X2RzbWTz8bLOZrDwpxy38ymLH+VTlUm2382+P6U55jd+90e649n0x1C6NpLtH0z3it/XxsDtcuqLdGV8O48hN4yWCyLyouaQz2EwYfuiP7Y/+mP7Ywu2LwxD0bLJhdQhmoceti/6Y/ujP7Y/tmD7AmvZWCwWi6UAWMvGYrFYLHnHKhuLxWKx5B2rbCwWi8WSd6yy8QkR2U9EfiEiD4rIxcWWp5iIyMkicq+I/FFEji22PMVGRPYQkftE5MFiy1IMRGS4iNzvXBNnFlueYjNUrwerbAARmSMia0XkXwP2Hy8ib4rIOyJydbo6VPV1Vb0IOBUoWzdHn/riYVX9FnAOcFoexc07PvXHUlU9P7+SFhaP/XIK8KBzTfxHwYUtAF76YzBeD26wysYwlwFRrUUkANwJnADsD5whIvuLyAQR+fOAz07Od/4Dky67ubDi+8pcfOgLh2uc75Uzc/GvPwYTc3HZL8AoYKVTrK+AMhaSubjvjyGJDcQJqOozIjJuwO7DgXdUdSmAiPwOOElVZwEnpqjnEeAREfkL8EAeRc4bfvSFiAhwG/CYqr6cZ5Hzil/XxmDDS78A7RiF08ogfcH12B+vFVi8kmBQ/vA+sTtb3sbA3DC7pyiLiIRE5A4RuQd4NN/CFRhPfQFcBhwDfE1ELsqnYEXC67XRICK/ACaLSOpMeOVPqn55CJgmIneT/5BHpUTS/hhC10M/rGWTGkmyL+UKWFUNA+F8CVNkvPbFHcAd+ROn6Hjtj/XAYFS6A0naL6q6ETi30MKUAKn6Y6hcD/2wlk1q2oHRCdujgFVFkqXY2L7oj+2P5Nh+6Y/tjwSssknNC8DeIjJeRKqA04FHiixTsbB90R/bH8mx/dIf2x8JWGUDiMhvgQiwr4i0i8j5qtoLXAr8DXgdWKCqrxZTzkJg+6I/tj+SY/ulP7Y/MmMDcVosFosl71jLxmKxWCx5xyobi8ViseQdq2wsFovFknessrFYLBZL3rHKxmKxWCx5xyobi8ViseQdq2wslgIgIpeLyOsiMr/YslgsxcCus7FYCoCIvAGcoKrLEvZVOgv/LJZBj7VsLJY840T43QOTfiIqIrNF5HFgnoiMFJGFIvKC8znK+U6DiDwuIotF5B4ReU9EdizqiVgsOWAtG4ulAIjIckwG10uBrwBTVPVTEXkAuEtVnxORMcDfVHU/EbkD+FBVbxKRLwN/Bkaq6ofFOgeLJRdsigGLpfA8oqqfOv8fA+xv8s0BsK2IjACOxqRTRlX/IiIfFV5Mi8U/rLKxWArPxoT/K4BggvIBwFE+dtjBMmiwczYWS3F5HDO0BoCITHL+fQY409l3ArB9wSWzWHzEKhuLpbhcDhwqIm0i8hpbMjjeCBwtIi8DxwIriiWgxeIH1kHAYikD4g4G1kHAUq5Yy8ZisVgsecdaNhaLxWLJO9aysVgsFkvescrGYrFYLHnHKhuLxWKx5B2rbCwWi8WSd6yysVgsFkvescrGYrFYLHnn/wMzg9s2/QWr/AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "res = ImpOLS\n",
    "rho = 1e12 * np.abs(res.impedance)**2 / freq[:, None, None]\n",
    "rho_err = 1e12 * np.abs(res.error)**2 / freq[:, None, None]\n",
    "phi = np.angle(res.impedance, deg=True)\n",
    "rad_err = np.arcsin(res.error/abs(res.impedance))\n",
    "rad_err[np.isnan(rad_err)] = np.pi\n",
    "phi_err = np.rad2deg(rad_err)\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = plt.subplot(2, 1, 1)\n",
    "ax.set_xscale(\"log\", nonpositive='clip')\n",
    "ax.set_yscale(\"log\", nonpositive='clip')\n",
    "ax.errorbar(freq, rho[:,0,0], yerr=rho_err[:,0,0], fmt='k.', label=r'$\\rho_{xx}$')\n",
    "ax.errorbar(freq, rho[:,1,1], yerr=rho_err[:,1,1], fmt='g.', label=r'$\\rho_{yy}$')\n",
    "ax.errorbar(freq, rho[:,0,1], yerr=rho_err[:,0,1], fmt='r.', label=r'$\\rho_{xy}$')\n",
    "ax.errorbar(freq, rho[:,1,0], yerr=rho_err[:,1,0], fmt='b.', label=r'$\\rho_{yx}$')\n",
    "plt.xlabel('freq')\n",
    "plt.ylabel(r'apparent resistivity  $\\rho$ ($\\Omega.m$)');\n",
    "plt.legend()\n",
    "\n",
    "plt.title('Site 002 results in 2-stage RR OLS\\n  configuration with sites 99 and 100 as remote references')\n",
    "ax = plt.subplot(2, 1, 2)\n",
    "ax.set_xscale(\"log\", nonpositive='clip')\n",
    "ax.errorbar(freq, phi[:,0,0], yerr=phi_err[:,0,0], fmt='k.', label=r'$\\phi_{xx}$')\n",
    "ax.errorbar(freq, phi[:,1,1], yerr=phi_err[:,1,1], fmt='g.', label=r'$\\phi_{yy}$')\n",
    "ax.errorbar(freq, phi[:,0,1], yerr=phi_err[:,0,1], fmt='r.', label=r'$\\phi_{xy}$')\n",
    "ax.errorbar(freq, phi[:,1,0], yerr=phi_err[:,1,0], fmt='b.', label=r'$\\phi_{yx}$')\n",
    "plt.xlabel('freq')\n",
    "plt.ylabel(r'phase $\\phi$ (degrees)');\n",
    "plt.legend()\n",
    "plt.ylim(-180, 180)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now computing 2-stage M-Estimator Transfer Function for site004 with sites 100 and 99 as remote references "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SignalSet: 9 channels, 1 run\n",
      "tags: {'site004_Ex': (0,), 'site004_Ey': (1,), \n",
      "       'site004_Hx': (2,), 'site004_Hy': (3,), \n",
      "       'site004_Hz': (4,), 'site099_Hx': (5,), \n",
      "       'site099_Hy': (6,), 'site100_Hx': (7,), \n",
      "       'site100_Hy': (8,), 'B': (2, 3), \n",
      "       'E': (0, 1), 'Bremote': (8, 6, 5, 7)}\n",
      "----------  -------------------  -------------------\n",
      "  sampling                start                 stop\n",
      "       128  2016-04-28 19:00:00  2016-04-29 07:00:00\n",
      "----------  -------------------  -------------------\n",
      "starting frequency 0.001\n",
      "starting frequency 0.00139742\n",
      "failed to converge (maxit=100). while processing step 2 (weighting=1).\n",
      "starting frequency 0.00195279\n",
      "starting frequency 0.00272887\n",
      "starting frequency 0.00381338\n",
      "starting frequency 0.00532889\n",
      "starting frequency 0.00744671\n",
      "starting frequency 0.0104062\n",
      "starting frequency 0.0145418\n",
      "starting frequency 0.0203211\n",
      "starting frequency 0.0283971\n",
      "starting frequency 0.0396827\n",
      "starting frequency 0.0554535\n",
      "starting frequency 0.0774919\n",
      "starting frequency 0.108289\n",
      "starting frequency 0.151325\n",
      "starting frequency 0.211465\n",
      "starting frequency 0.295506\n",
      "starting frequency 0.412946\n",
      "starting frequency 0.57706\n",
      "starting frequency 0.806396\n",
      "starting frequency 1.12688\n",
      "starting frequency 1.57472\n",
      "starting frequency 2.20055\n",
      "starting frequency 3.07509\n",
      "starting frequency 4.2972\n",
      "starting frequency 6.005\n",
      "starting frequency 8.39151\n",
      "starting frequency 11.7265\n",
      "starting frequency 16.3868\n",
      "starting frequency 22.8993\n",
      "starting frequency 32\n",
      "(32, 2, 2)\n"
     ]
    }
   ],
   "source": [
    "from razorback.weights import mest_weights\n",
    "from razorback.prefilters import cod_filter\n",
    "\n",
    "sig = prepare_signalset(inv, 'site004', ['site100', 'site099'])\n",
    "print(sig)\n",
    "ImpME = rb.utils.impedance(\n",
    "    sig, freq,\n",
    "    weights= mest_weights,\n",
    "    remote='Bremote', # including the remotes references in the computation,\n",
    "    prefilter=cod_filter(0.0), # no coherency prefilter...\n",
    "    fourier_opts=dict( Nper= 8,  overlap= 0.71) # fourier options with 8 periods by window, and 71% of overlap\n",
    ")\n",
    "print(ImpME.impedance.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-17-a0bace2a0562>:5: RuntimeWarning: invalid value encountered in arcsin\n",
      "  rad_err = np.arcsin(res.error/abs(res.impedance))\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(-180.0, 180.0)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZsAAAEpCAYAAABFmo+GAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABfmklEQVR4nO2deZwU1bX4v2dmmEVmGHUAFxbBBZ8kIBiM9gvRNuMaTTRqXBJ/uAY1LjExJu67oHm+xBCXiA8ekmASE4kxLlEzsVVC8xRlHCLiBgKjyBZph222Pr8/bjX0DL1Ud1evc7+fT3+6q+r2rVO3q++pc++554iqYrFYLBZLNinLtwAWi8ViKX2ssrFYLBZL1rHKxmKxWCxZxyobi8VisWQdq2wsFovFknWssrFYLBZL1rHKxhIXEfmuiLyQbzkKFREZISIqIhUe1PWciJzrhVwWSyFilU0fR0Qmish8EQmJyL9F5J8iciiAqs5R1WOjyqqI7J/BuUaIyEsiskVElorI0b2Of0dEVojIZhF5UkR2jzp2r4i8LyJtzncnpStHthCRgIhclM53VfUEVX00jXNWicgMp93aRGSRiJyQjgwicquI/Dad72aKiJwnIt0isklEPheRt0TkpKjjEcW+yXl9JCLXJqlTRWRN9MOAiFSIyFoR0ah9ARHZFlX3JhH5a3autO9ilU0fRkQGAE8DvwJ2B4YAtwHtWTrl74BFQANwA/AnERnkyPIF4GHg/wF7AFuAB6O+uxn4BlAPnAv8UkT+081JvbA8CpgKYBVwJKZtbgIeF5ER+RQqTYKqWgvsivntfy8iu/Yqs6tT5nTgJhE5JkmdG4Fo5ft14LMY5S5X1dqo1zfSuQBLAlTVvvroC5gAbExw/DxgnvP5FUAxnf4m4Exn/0lAM+ZPPR8YG6euURglVhe171XgEufzFOCxqGP7AR3R5XvV9xRwdZxjfqAV+CnwKfAbzIPVtcCHwAbgcWB3p3w18Ftn/0bgdWAP59hHwNFRdd8K/Nb5PMJpkwrgLqAb2Oa0z/2AAL8A1gIhoAX4YhyZA8BF0e0O3IvpGJcDJ6Twu7YApyU4/lPgY6ANeBdoBI532rvTkf8tp+z5wDtO2WXAxb3q+gmwGvgEuMhpj/2dY1XONawE1gC/BmqS3WvO9i5OXYf2buuoMq8B1yS4TgVuBP4Yte9PmAcdjdX29pW9l7Vs+jbvAd0i8qiInCAiu8UrqKpHOB8PVvPk9wcROQSYCVyMsVYeBp4SkaoYVXwBWKaqbVH73nL2R46/FXW+DzGd36jeFYlIDXAo8HaCa9sTY63tA0wGrgROwVgAe2M68QecsudirIJhznVcAmxNUPdOqOoNGOUZeUK+HDgWOMK5hl2BMzEKzQ2HYRTBQOBnwAwRkWRfEpE9nPPFbBsRORC4HNOJ1wHHAR+p6t8wCv8PjvwHO19Zi3mgGIBRPL9wfndE5HjgR8DRwP6Yto3mHkeWcc7xIcDNLq6h3DlXJ7AiTpnDgS8CHySp7kngCBHZ1bGSvgr8JZkMFu+xyqYPo6qfAxMxT4CPAOtE5Cmnw3LD94CHVfX/VLVbzZxDO3B4jLK1mKf7aEJAncvj0fwao5ieTyBbGLhFVdtVdStGId6gqq2q2o6xUE53htg6MUpmf+c63nDaJlM6Hfn/AxBVfUdVV7v87gpVfURVu4FHgb0ww4txEZF+wBzgUVVdGqdYN8biGC0i/VT1I0exx0RVn1HVD9XwMvACpsMGOAP4X1V9W1W3YIZgI7II5v74oar+23nImAKcleASDheRjRjr8F7gHFVd26vMehHZCgQxQ21PJqgPp66/YhT9WRiLeFuMctNEZGPU644k9VpSxCqbPo7TAZ6nqkMxT4p7A/e5/Po+wNXRf1KMdbB3jLKbME/H0QzADM+4OQ6AiPyXI+cZqpooiuw6VY3uVPYB/hwl5zuYjncPzDDb85g5gk9E5GdOx50RqvoPzHDaA8AaEZnuzJO54dOoerY4H2vjFRaRMsx1dGAsl8j+56Imvb+rqh8AV2GU7VoR+b2IxPq9It8/QUQWOM4jGzFzHgOdw3tj5osiRH8ehBkKeyOqzf/m7I/HAlXdFdgNoxS+GqPMQEw7/BgzXOrmd5oNTHJes+OUuVJVd4163eSiXksKWGVj2Y7zNDwL05m7YRVwV68/6S6q+rsYZd8G9hWRaEvlYHYM97ztbAMgIvtinsDfi9p3G2ay91gXlkdvRbQKM+8RLWu1qn6sqp2qepuqjgb+EzNsFPF224zpNCPsmcI5UdVpqvolzDDhKOCaJHKnjGNFzMAoztNUtTPq/CfojknvOc6+x1R1IkYBK2a4ayf5neHQJzBWxh6OIngWMxcFZq5maNRXhkV9Xo8ZivxCVHvXq5ncT4iqbgK+D/w/ERkf43i3qv43xkL5frL6MMObEctwnovylixglU0fRkT+Q0SuFpGhzvYw4GxgQZyvrAH2jdp+BLhERA4TQ38RObGXQgFAVd/DOBLcIiLVIvItYCymMwMz/PMNEfmqiPQHbgfmRuZ4ROQ64DvAMarqdt4jml8Dd4nIPk59g0TkZOfzUSIyxpkr+Bwz/NXtfK8ZOEtE+onIBIwXVDx6tI+IHOq0TT+M0toWVa+XPAQcBHzDGTKMi4gcKCJfcxTJNoxCiMi0BhjhWEkAlRiFvw7oEuNSfWxUdY8D54vIQSKyC1HzMaoaxtwfvxCRwc65h4jIcW4uyPmN/4fEczx3Az8RkeokdSnGk/GbSaxhSxaxyqZv04aZiP4/EdmMUTL/Aq6OU/5W4FFnWOQMVV2IGZe/HzPh/gHGqygeZ2E84D7DdBSnq+o6AFV9GzMxPwczKV1Hz6fWKcBw4P2oYaHrU7jWX2KGZl4QkTbnWg9zju2J8VL6HDO89jLGOw2MK/F+jsy3AY8lOcfpIvKZiEzDDAM+4nx3BcY54N4UZE6KozwvxkzCfxo9ZBbnK1WYtl+PGaobDETa8Y/O+wYRedNR9FdilMpnGGX/VKQiVX0OmAa8hPntg86hiOv8T539C0Tkc+DvwIEpXN59wNdFZGyc4884cn0vWUXOvFIih5L7pec6mzdSkNPiAslE0TtPoNucSUyLxdKHEZGDMA8rVaralW95LIVFSpaNiJSJWeX9jIisBZYCq0XkbRH5LxE5IDtiWiyWQkREviUilWLc5u8B/moVjSUWqQ6jvYQZUrgO2FNVh6nqYIzXyALgbhE5x2MZLRZL4XIxZk7nQ8zcz6X5FcdSqKQ0jOb45XdmWsZisVgsfYu052ysUrFYLBaLW9LyRhOR/8EsUlslIv8nIo+IyBUey1YUiMgeIvKKmIi7/y0i1zvtUxCIyK9FJG8L1JKdXzyKNJzv68w3YqIm2zUkRYqI3Cki60Xk0+Sli5N0XZ+/ilnkNQw4Ffgz0N8zqYqLyRg30gGqerWqTlHVtMLMZ0qsDkdVL1HVvIXeiD6/iPhFpLWYz+OsKfmHmJQMHzjrhaKPX+Ts3yQif5MEq/NziZjoBe+KSFhEzotx/Ici8qlzXTMlKr6diOwuIn8Wk/phhYh8J6fC5xivHoBSON8wzHKD0aqaaNFwUZOuslmACSmBswL7WVW92zuxiop9gCXZXiwmpR0mvyhwfoO/YNIy7I550PitiIxyjh+JWQ90snN8OSatQiHwFmbd0pu9DzgLLa/FRH8egVmYeltUkQcwYXD2AL4LPCQmJUTOcBbcFiUu/rv7ABtixIFzU7dELcItbDSNUNHAVzCL336MsXLq06kn1y9M0Mn5mDDyq4DznP31mJhJ6zCL724Eypxj5xEn3DsmtEsn5o+4CRP99lacEPROmUnsWNB3E1Eh653v3xlV1g+0Rm1/hFkY14JZKFfBjjD5bcAS4FtO2YPYsUJ9E07qgBjn+B5mod2/MQv09o46ppiFle+zIyqyxGjHaszK84HO9o1AF8a6A7gTuC/6/BjLdysmQOYm57W3016PO+3fhglbMyHO7xc3ZL+L86SVYqDX+b/o1CdR+14A7nA+3ws8EHVsb6dN94tzPXHD97MjTcLVzvWuBs6POt7g/H6fY0Lt30FUiP4E/4F5OPd91L7HgClR243Ap87n/jjRt6OO/wa4O079+wH/cNpyPWaR7q5Rx3dKbxCnnlmYyAjPYqIvHO205xOY/+lyTDyzSPlbMYtSf+vUvRgTHug6p/1WYcIcRf82T2H+Bx8A33P2x0u1UI8JCbTakf9OoDyO7LdiFgn/1vl9Lor3fee6ou/XWU4dh7Ojr3oL8EfVH8CktPin8939McFeX3Su511M7MDotnwAswi2Dfg/ou5JTCilyHfXANc7+zP+z/RolzQ77eWYJ59rMU9ubwMfplNXrl6Y1edtmHAs/TB/1nHOsdmYJ9Y6zJPde8CFzrHznBvve87NcSkmd0fEuWIWPTvzW9mR72S0cwNNxIT+uNepKxVl04yJOVXj7Ps2OzrPMzF/xL2iZJ3X67q3nwP4GqYDOASzkvxXwCtRZRXz1L6r017rgOPjtOcrODlTMB3uh+xQwq+wQwlGn7/H9UW11zZMgMdyYComIGOscx4HvOHIJxgFu5fL81yFsciHOtf+MPA759jFmMjAuzgyfAlHcfaqYww7K5sXgT87n/8beDDq2BCnTU+Ocz0nYjpnwYTn3wIcEnUNXZiwPf2c9tkC7OYc/z3mz98fowQ/7v3bxzlnLGXzFk5+Imd7oCN3AzAe2Nqr/I8x62li1b8/cIzTxoOce+E+59iBmE5/b2d7BPEV8SzMA8VXMPf6Ls5vfzPmv7QvRkEf1+s+Og7zUDYb00/d4LTf94DlUfW/jIkaXY2JvrAOR/HR64HR2fekc8/0x0RdeI1euX163dOdmJQWZUBNou+z8/9+CKYT/7rz/WOc7UHO8QAmR9AXnGutd9r1fGf7EMz//AtRbflv4MvO8TnA751jdRgFeLXTFnXAYV79Z3q0S6qdtnOiV2Psq0qnrly9ME84f46xvxxjNYyO2ncxEHA+nwd8EHUsktRpz6gfMp6yuTny40R9t4PUlM0FSa6rGaczI7mymQH8LOpYLeZPMcLZVmBi1PHHgWvjnPcOTKiSCkzYkx9gwqD0tnqiz9/j+qLa6+9R26Pp1blFHfsa5kHgcBzLM851xjrPO0Q9RWMCM3Y68l9AgsRvUd/ph+ngfuJ8Ptb5PZ93jjdi/uRjMR3Mw5gn1rNd3qNPAj+Iuoat9EwWtta59nJH9v+IOjal928f5xyxlM2HRD1UONemGGXwVRwrJ+r493D+Hy7OdwqwyPm8v3MNRwP9knxvFjA7avswYGWM//T/Rt1HL0Yd+wbmwaDc2a5zrmlXzMNbNz0T+U1lh1VxKz1HJ/bA9BE1UfvOBl6KI/ut9HyIS/j93vcrxvr7Ta86nwfOdT4HgNujjp1Jrz7ZufduiWrL/4k69nVgaZQci+JcR8b/mehXumN9zSLyg+gdanKEFDLDMH+q3gzEPClFJ2lagXm6iJBSuPcoeoRgd76bahDJ6LDtiMgkEWmWHWHbv8iOkO9u5Nl+nWqi624gzrVinqTjXefLmD/JIZghixcxT+eHY5TzepcyxTpndaxxbs0sZP8+ZJhiQI2r/ykYi+RTzNPg45jhLlS1CbgFM9SzAvOw0BY53htJHL4fzDh+9Gr8yO8xiB3poCPETDLmkt7pHSKf22IcixxvIwYiMlhM2oKPnXhov8W5Jk0xvQE9r28fYG/pmc7ienrm+FkT9XkrsF53hNKKBCitxfwPIjl2IvT+z0ezD0YBr44698MYC8Wt7Kl8fx/g272udSKms49X/2G9yn+XnhHK4/2v4/WLkXo9S8uRrrLZAxPt9xMReVpE7hKRb6dZV65YhRmy6M16jLbeJ2rfcMywRKb0CMEuJsNkQ9RxN+HrNer7+2ACO14ONKgJ+f4vdoR8152+3ZNPiLpOMbHtGkjvWudjhkW+Bbysqksw7XYiRhHFIpl8SVF3IftjnSfdFAO9z9+iqkeqaoOqHocZznkt6vgDqnqAmsgaT2CUwr961yPJw/cnYh1miC06pP9wF9+LR4/0Ds7nNWoiL78HVEjPUFQHEz9L6lRM+49V1QHAOURdk8ZPbxCL6N9xFWYYLPr3q1PVr7u7xB58AuwuPaOTR//ne98/qzCWycCocw9Q1UROEr1lT+X7qzCWTfS19teeTli963+5V/laVXUTzSFevxg5lvF/JkJaykZVz1DVg4CRmKGi9zDjgYXMHOBoETlDRCpEpEFExjlPPo9jws/XOR36j9gR9TcT/oQJm/+fIlKJmeeK7kyaMVFtdxeRPTFPfYnoj7nJ1gGIyPn0zD2zBhjqnCsWj2FCwo9zOrspwP+p6kcpXRXbrbQ3gMvYoVzmY4Yg4ymbNUCDiNSnej5IKWR/rPOkm2KgtwxjxaRI2EVEfox52pzlHKsWkS86HkLDgenAL1X1sxhVJQvfHxfnnp0L3OrIMRqT2jouYuKXVWPuv36OrJH//2zgQhEZLSbG2Y2Ra1LVzc65bheTQuIrGG+738Q5VR2Og4qIDCHqYUASpzdIxmvA5yLyUxGpEZFyp60Pdfn97ajqKsy9OtVph7HAhZg+AnqlWlCTXfUF4L9FZICYGJH7Od6Hbs6X6vd/i+k3jnOus1qMO//QOOWfBkaJyP8Tkwqjn/NfOciFeE8De4rIVSJS5fSBkWjonvxnImTkMqcm5e6bqvqoqnqeFMpLVHUlZpjiasxkWTM7nuauwHReyzBj2o8BMz0459tO3b/HWDltmDHryJDjbzCTsx9hbsY/JKlvCWYSOoj5Q4zBeKRE+AfmifNTEdlpGMsZ5rkJ80S9GvNEkyhNbzJexgwPvBa1XYeZFI4l/1KMQ8kyxzRPdQ2Kq5D9cc6TboqB3vw/TNutxczRHBM1hFyNuXc2YdokiGnvndAk4ftdcDlmKORTjGL43yTlX8B07v+JUYJbgSMcWf4G/AwT+3CF87ol6rvfx8xBrcW066UaP1z/bZih1RDG+2lu1LFE6Q0S4ijYb2Am85c7dfwPZnI8Hc7GzEl9glkneIuqvugc65Fqwfk8CfOAsATze/2JnsNayXD9fUcZnoxpm3UYC+Ma4vTXzr10LOa//Ammbe/BtHdCnO8eg2nbTzGeqEc5h736zwAZphiwpIaI1GLcBA9Q1eV5FsdisVhyhmeLgZxhIEsvROQbzlBHf8xT+GKMJWOxWCx9Bi9Xns7wsK5S4mSMafsJcABwllpz0mKx9DEyifpcEe2WKSL7qGom7pcWi8ViKVFStmxE5Hsi8i6wypl8/YeIHI6ZZLNYLBaLZSdSCu4oIj8FJgBHquqnzr5jMR5Cmfj5p4WInIJZ1zEYE5PqhVzLYLFYLJbkpJqpcykwRnslTRORPTAxoU7LWCCRmZgFQmtV9YtR+4/HuOKVY0Iv3B11bDfgXlW9MFn9AwcO1BEjRmQqpsVisfQp3njjjfWqOijd76cctr63onH2rRGRX6UrRC9mYUKSzI7scBYOPYDxB28FXheRp5x1J2AWoT3gpvIRI0awcOFCj0S1WCyWvoGIZDQnn+qczYcicmIMIW4HmjIRJIKqvoJZdBnNlzHxtpapagdmkeTJzkrte4DnVHWnPB2eEgzC1Knm3WKxWCwpkapl833gCTGZ/t7CrGD+OmY1/rueStaTIfQMPNeKWcl6BSaCbL2I7K+qv471ZRGZjEl0xfDhaUwtBYPQ2AgdHVBZCU1N4POlXo/FYrH0UVKybBzX5kMxa2o2Y9aOnKOqk4BHvRdvO7GCE2okKKOalMAxFY1TcLqqTlDVCYMGpTHkGAgYRdPdbd4DgdTrsFgsFo9xO+ASnL6YqccFCE5fnBvBYpCqN5o4CxL/5ry2o6r39CrjJa30jHA7FKPocoPfbyyaiGXj92deZzBolJbfb60ki6UEcPuXdl1u+mICT2zAf1oDvsljYtbTeFQ3HR1CZaXS9FJ5zPqC0xfTePF+dHAQlS900MTimPVlm1SH0V4SkSeAvziBLQETURaTb+FcTDC/WZ5JaHgdOEBERmLCgJ+FCVyYG3w+M3TmlXJwOyxnFVJREQwGCQQC+P1+fPb3Khnc/A1dd/weKojA7BV0tA+hm3I62jsJzG7F59tnp7oCT2ygg4PopoIOlMATG/BNTr0dMiVVZXM8JkPb75yOfyMm0m05JqrsL1S1OROBROR3mKRcA0WkFRONdYaIXI5J1lMOzEwQdTY7+HzedfixhuV6123niYqKYDBIY2MjHR0dVFZW0tTUZBVOgeOlEnHd8XuoIPy8TCWn04FSSSd+XiZWShn/aQ1UvtCxo9xpDTuVyQUpKRtV3YbJ2/2gk1NkICaF70avBFLVs+PsfxaTXKr4cTMs50YhRbAWUN4JBAJs3WqSQXZ0dBAIBKyyySPJ/hJeKxHXHb+HCsI36QCaZn6dQOdX8Pf7J75JU2O2hW/yGJpIPCSXE9Rl/uhSeX3pS1/SgmD+fNUpU8x7vOM1Narl5eY903KWrDJ//nytqanR8vJyramp0fn2d8gKyf42kTI1VV1aLt1aU9UVs+yUSz7ScjoVVMvp0CmXfBS7rkse1Ro2azkdWsNmnX/Jo3FPOr/ySJ0i1+v8yiMT/l9dlVPV+Q+36JRjX9L5D7ckvthkDeIRwELNoO/Ne+ef61fBKBs3uLmRpkwxigbM+5QpmdVnSZv58+frlClTrKJJA6+UiKo7ReK5EnF7EamUKzCssillZeOGVCygsjLzk1sLKGXy0T/0FeXl9hb21BrJhhIpcTJVNimHq4ng5Ocu06g0A5Y84NZTLhAAcZYreTUH1EfmioJBmDgRwmGoqcmNr0YwGOQrX/kxqkdQWXkdgcDUopwDcnOLBALQ0a50h4WOdiUQkJhl3c53uJrL8PnwBabiCwTAPzXxD+qlc1BfJh0Nhcl/vg7jhvwWcFEmGi+Xr5KzbNzi9RxQH5ormjLFGIRuRiq94pJLHlXYrNCpsFkviTfMU8C4Hfaa/3BLT0sk3hyFtUbyCnmybK7GRH/+VET2AqaIyFBVvdUD/WfJBqlYQG684NyWKwErye83Fo2Xa3qTcyRQiXEYVWc7+3i14BDce3L5NjxNU9kzBMJfxV/2Kr4NJwIx6rTWSHGTjoYCFmOG0CLbZUBLJlovV69CsWzcjMfnZczeS8vG7TxRPq0kl0/A48ZdqiNHPpwz0ebPV62q6lKRLq1KYBV4fU5XP/12S6QzoSWSyiR8X7GSixnyZNk8BPxRRH6qqh9gEqdt8UD39QmCwSATJ04kHA5TU1MTcwFg3hYJurWA3JRzO0+UrzVFKSycra9fQn39Enw5Wnrt88FLL5V7auwlazq3cyduV6S7XQfieYQOS0GSlrJR1QdF5FPgf0RkLDAAmCMi3waaVfV9L4UsNQKBAOJ0wvEWAOZzkWAQCGDCOCQ8Y7KhCrcx5dyW8zqqQiAAThsnVXIe4neuL5AkoKuXI0Fums7fsJjK8H500I/KcCf+hg+JNZzlekW6HfayRJG2N5qqzgXmikgFMBoYD/wncCnwNW/EKxy8jHvl9/uprKzcbrX4Y3Sufr+fmpqahGWyIZunFpWXVhKkZgG5we8nWHmkefIu/ye+DCdj8jXt5Ob3d9N0budOUlqRbpWIJUImY3DF+Epnzmb+/PlaVlamgGerw72as/F65fqUKVMUMyOt5eXlOiUXrlduSWVs38VczPz5qhXSpUJibylV1XHjxunIkSPjtm8qoh155JF65JFHxi+QAvPnz9fKyiNV5HqtrDwyM/ns3IklAeRrnU0p0dnZSWtrK9u2bYt5XER45plnemy/8847GZ1z11135ZRTTgGIW5ebMiLC3LnPAFWEw9vYtGkTnZ2d9OvXb6eybp6AU7GovMSVVeDzEbzv/3Y8UfviPFG7HG4LBKBLywHo6ErkUBekpaWFcDhMY2NjTGsvFaNr4cLpdHXVEgwmeeh30SizZ79PR8ezgLGUZ8/+U8zf1pXxaOdOLFnEKhugtbWVuro6RowYsX0uJZpNmzbx3nvvEQ6HKSsrY9SoUdTW1mZ0zk2boK0N6uogk6rWrdvCihVVGIfAbnbbbQ2tra2MHDmyRzk3TgkAPp+PpqYmD4flXEbWdZlxofGqMabcq9A0Jk6dLnt+ty7NgUCAcDgMxJ8/S2XaafPmfQHhqKO6eSlOAEj381PuXaRdjWjZYS9LlkgpU2epsm3bNhoaGmIqGoDa2lpGjRrFkCFDPFM0770HH39s3jdtil9u9er4xwG6unbB/IwClFNTs3dMC810mF8GrqW9/ZAkk9M+4DqSuAckzRIY6S9vusm8xysXmadPlgjVdcLUSM9fXp6w5488yN9xR2JfA7/fz3+KcC0wsbw8prXn80HTfYu5ozFA032L49Y1e/YK51M57e3hqO0YF+uiUSZN2oeqKkGkm6qqMiZN2nkdi8VSCFjLxiGeoolQW1vbQ8m49SiKRVubCX0C5r2tbWfrJqKQwmEoK4NRo2JbQHV1UFYmTjlhwAAIhXYu19BwEvADoJJwuIOGhg9jypaSlZGknNuhJbcWhuuEqW6H23D3IO8DXlY1uclFKI9VKBjEd+lEfOEwvFoDY+Jpr5eB0zFWSKezvXPIFbeNkpKLdAEvnLWUPlbZ5AGjIHYokrq6ncu4UUhg9o0alXxIbsOGMZgOTigrK3e2d8bLAAJulYPXTmuuh9vcEghQjrEd6YozuRMI7PjBEjTcpEkH8PDDx2BingWZ5MHaE1cjXzYZnyXPpKVsnKyZc1T1M4/l6RO4URBuFFJ0fclG9syDsjh9jWRsPfj9EA5vAyqoqBD8/p2f91OZb3Y7VeCmnOfLZ/x+2svKqAiHqYjXKK4tER8HH7yFUOgPzJkzJ/GcmNtGcR3t0kO3cYslRdK1bPYEXheRN4GZwPOOa1zJ0ntCPxQKEQqFCAaDaU2iJ1MQbi0Wt3htPZilnz8GjkA1CEwl1hzP978/nlAohN8/J+bxbOB5LDOfjx+NHcu4UIjJc+bEbpQUNGt9fT319fXeLNJ1a7G4HoO0WLJEuj7TmFGF44DfAx8AU4D9MvHDzsUr1jqbJUuWJPQvb2tTfeMN1ddfN+9//7v3624iPPDAA/qFL3xBhw8frtOmTUurjmTX4wVu1uNkY32SW9wG/XW75sXLtTFe1mWT51lyBflaZ6Oq6oSs+RToAnYD/iQiL6rqTzJRgIVG7/mTpqbkrrDp8MQTT/Diiy+yaNEi1q9fz5gxY7j00kupqCi8qTW/309ZWRnhcDjuehw3LsPZopA9eEeHQowLhUi+0MYFqVgshdwolpInLddnEblSRN4Afgb8E5Nu4FLgS8BpHspXEETmT8C8NzaajhbwdOHjtGnTuOeee+jXrx977bUX/fr1295ZFxo+n4+xY8cycuTIuGt2IgoJkrfT+PHj2XfffQnG848uFYJBft7SwoXLlyf2B3eLW/9tiyXPpLvOZiBwqqoep6p/VNVOAFUNAyd5Jl2BEJk/GTLEvDc2Ju9oU6Wzs5OWlhZGjRoFwOrVqxk4cCChUIjzzz+f1tZWLrjgAtatW9dju7OzM+Nzp8uiRYtYtmxZ3Ot3o5Bgxwr95cuX09jYmFDhFL1SCgSoCoeN+3TCxUIp4PPBddclVTTBYJCpU6d60nZe1mXpG6Q7PlOlqj1Wo4nIPar6U1XNLI5LgdJ7Qt/TSV5gyZIlhEIhli1bxogRI7juuuu44oorGDRoEMOHD+fqq69mxowZ1NbW9tiOFZamkHDTTm6H29yEjYmQyTqorOLGsy0LuI0gkeu6LH2HdC2bY2LsOyETQfo6ixYt4rvf/S5nn302Y8eOZfjw4UyePJlNmzaxbNkyKioqqK2t3Wm70AkEAkk7fLfDbbGUUtHheLbNHDkyp8NeXrZdrBQZFksyUrJsRORS4PvAviLSEnWoDjN302fw+g/W3NzMSSedxJlnnrl9X1dXF1deeSV33nknjz/+OH//+9957LHHtm9H4pcVO5HhtlAolHDtiRunhFTJ1IU9HZbU17Okvp7JObQGvGw7NykyLJbepDqM9hjwHGZRxbVR+9tU9d+eSdUHaW5u5tJLL+2xr6KigpkzZwJwzTXXAHD00Uf32C4VFi1alLSMz+dj3rx5nubucTssV+x42XZeB2u19A1SUjaqGgJCwNnZEafvYoci3OHz+Tzr3PLpmp0PvGw7L+uy9A1SmrMRkXnOe5uIfB71ahORz7MjosWSOqFQiJUrVyb0lkrFNbuvYL3MLNkiVctmovOeIFJX7hCR/sCDQAcQUNU5eRbJUgC4HR5zO1fUV7BeZpZsku6izh+KyBCvhXHqnikia0XkX732Hy8i74rIByISmS86FfiTqn4P+GY25LEUH6l4XtXX1zN8+HDbqVIi3n6WgiVd1+cBwAsi8qqIXCYie3go0yzg+OgdIlIOPIBxrx4NnC0io4GhwCqnWLeHMliKGDs8lh623SzZJC1lo6q3qeoXgMuAvYGXReTvXgikqq8AvT3bvgx8oKrLVLUDE/zzZKAVo3Agx1lH/bP8+Gf5c3lKi0vcRi6w9CTisTZlyhTbbhbPyTTC41pMIM4NwODMxYnLEHZYMGCUzGHANOB+ETkR+Gu8L4vIZGAywPDhw7MopqVQ8DrCg5eEBoQI7RoiuCqIb1hhyWe9zCzZIt05m0tFJAD8AxMn7XuqOtZLwXqfMsY+VdXNqnq+ql6ayDlAVaer6gRVnTBo0KAsimmxJCa4KkjzIc0s33c5jbMbCa6K7/UVXBVk6qtTE5axWIqFdIeehgM/UNXRqnqLqi7xUqgYtALDoraHAp9k+ZwJCbWHWBla6XlH8OCDD/LFL36RffbZh1/96lee1m3JP4GPApHs3HR0d5jtGARXBZn4vxO5/h/XJ1VKFksxkNY6G+AK4NUcrrN5HThAREaKSCVwFvBUFs+XkOCqIC1rWli+MfnTaSpE57NZsGABt912G11dXZ7UbSkM/CP8lGkZhKGyvBL/CH/McoGPAoTV8QxLoJQslmIhJWUTtc6mVlUHRL3qVHWAFwKJyO8wOYcPFJFWEblQVbuAy4HngXeAx1X1bS/Olw7Z6giKKZ+NJT18w3yMbR7LyOUjaZrUFHfOxj/CT01FDeVSnlApFTp2kaglQloOApF0Asn2pYOqxgyFo6rPAs9mWr8X+Ef4KZMywhr2rCOIl89mzpw57Lnnnhx//PFceOGFHHrooYwYMWL79gMPPEBNTU3G57fkjvrP66n/vD6hc4BvmI+mSU0EPgrgH+EvOEcCN9hFopZobIqBNPAN8zF2j7GM3DXx02kiNnVsYnXbajZ1bAJ65rMJh8Pb89kcccQRzJs3jxkzZnDmmWdy7LHH9ti2iqZwGP/wePb95b6eDav6hvm47qvXFaWiAbtI1NITm2IgTeqr6qmvSvx0Go9NHZtYun4pAGWbyhjVMKpHPpvNmzdz6qmnMnnyZMBERA6FQlx00UUxty35JzKPF9YwjbMb034IKSX8fj81NTU2FYEFsCkGXLOpYxNt7W3UVdVRW5lZ0rK29rbtn8Mapq29LWY+mwgVFRXcfPPNcbctO7PwSwvp6teVdC1L87hmV/UlKxdrHq+vKxubisASjU0x4IJYlkjgvEDa9dVV1VG2ycz5lEkZdVV1MfPZfLzuY2684UbO/O6ZDB48mFAoxA033MC5557L4MHZXENb3ARXBdlctxkgZ1ZGNubxSgG7SNQSIV0HgW8Df1PVNhG5ETgEuFNV3/RUugIhliWSiXVTW1nLqIZRPSyl3uPZmzo2saZrDZfdehllUsamjk3U19dz//33p33evsJ270DJnZURmccLbQsx59Q5fd6qsVh6k264mptU9Y8iMhE4DrgXeAgTQqbkiGWJZEptZW1ChdXW3rZ9WMYLBdeX2G5VKFRW5M7KWHRx8myjFktfJV1vtEiE5ROBh1T1L0ClNyIVHhFLZEjdEEY1jMpJp19XVUeZmJ/HKwXXV/AN89G/rT9V26rsRL3FUiCka9l8LCLTgaOBe0SkihxHXc41ySyRbJyv91CbxT0T3pgAgO9uq2j6IsFgMKljgpsy+SxXaqSrbL6NyTnzM1XdKCJ7Aj/2TiwL5F7BWSyFzvjx45NmVnWzmNTtgtNgMEhjY+N29+1clStF0o2NthZ4FHjTiYn2PvCEx7JZLJY+gt/vT7oOJ5Lue/ny5TQ2NsYNgeNmMWkgEEBEEpaJlOvo6KC7uzun5UqRdGOj1WUrNprFYskOxZ6ywG1Egshi0vLy8riLSf1+P5WVlQnL5LNcKZJp8rS+S+QmKaInk+CqYFHH2rKkTyRlQVjD1FTUFJzjRCgUIhQKEQwG4w4ruY1I4GYxqdsFp/kqB6U3t+PFOpubgPGU8DqbQsdNdINC72ws6TP+4fFJ1/cUcoSDyPBYOBymsbEx7jxGKh21m8Wkbhec5qNcKc7tpOtBdpOjaCYCx2Lmbx7yTiyLWyLRDT5u+5j3Nry3PbBnb2x+lNLEbW6lQk5ZkErATp/Px3XXXVf0HW8y3M7tFFMKB7vOJl1CIVi5Ejz+kVPN1BkrukEs8tXZ+Gf58c/y7lxeR1Yudtw+RERSFtxx1B0FZ9X6/X7KykxX1NfmMeLhZm4n4lF3/fXXJ3SYiJTNt1LKZJ3Nw5hUA31inU0PgkFoaYFwGBoboakJPHjSis7UuX79esaMGcOll15KRUX8nylWdIONbNypXDbyo7gZvvESG1l5Z1KJyeYb5ivI9vL5fIwdOzapS3Nfws2QYSyPukJ2t05X2ZyBWWdzr7POZi/gGu/EKnACAaNoADo6zLYHP960adN45JFHUsrUmcriTy87m3x0/IU875AvfMN8zDt/XtE7ftTX11NfX28VTTRDgYnOewwi1k8yh4lAIED7oHbCw8O0r2yPq5SyTbrKZivQHxP9+XagH8R4nC5V/H4oKzMKp7Jyh2daBsTL1BkKhfjJT37CHXfcwc0334zP52Po0KE9MnXW1ni7+NON11o+On4bWTk2hWqxWNInuCpI4+xGOro7qCyvjPkw59ZhomFcA+H/F4ZyCHeHaRjXkItL2Il0lc2DQBj4GkbZtGEWdR7qkVyFjc8HY8eaeZs5czyxaqIzdY4YMWJ7ps5BgwYxfPhwrr76ambMmMGaNWuYOXMmH3/8cVYydbr1WstHx28jKxcOkXm4TFJt9FXceg92dHfQrd2JH+aSWD8AG2o3mJ5eoKyszGzngXSVzWGqeoiILAJQ1c9EpG85CNTXm5dH5mi8TJ2bNm1i2bJlVFRUUFtbS21tbVYzdQY+CiA448AJbvJ8Dd94HVm5q7zLVZI1i8UL3A4/+0f4qSyv3G7ZxHqYc2P9ROqq6VeTsK5ckK6y6RSRckABRGQQxtKxpEmsTJ1dXV1ceeWV3HnnnTz++OPbzeVsZup0c5NHKPbhm3wkWbP0bdwOP7tx6Al8FGBr11ZP6soF6SqbacCfgcEichdwOnCjZ1IVAx5HDoiVqbOiooKZM2cCcM011xAKhbj88suzmqnT6xsz1B4itC1UkJZDKknWrAVk8QIvvQcjyxmK5cEwZWUjxtfuFeANoBEQ4BRVfcdj2foUbgLy5SpTp1c3ZqG7KrtNsmYtoJ3x+iGirwSk9HL4uVAsFrekrGxUVUXkSVX9ErA0CzJZSoRUPNbyMeEcSbLW1a+LpgvjK5B8pJkuZAr9IcJrvI4p6KWVUQgWi1vSHUZbICKHqurrnkpjKSmKwVW5oruCiu6KpMMVQM7TTBcq+VzvlOuHEhtT0DvSXfV/FBAUkQ9FpEVEFotIi5eCWYqfiKvyyF1HFvWf1KaZ7knkIQIo2IcIr4jlnWlJj3QtmxM8lcJSstRX1VNfVZ+0gy5kRwJwZwGVAm4sh2ysd/LaYvGqvlS8My2JSUvZqOoKrwWx9F362hxAKeD2IaLYKbZJ+EKmqJOnicgpmMjTg4EHVPWF/EpkSQcb86x0KeRIA25lK6ZJ+EImb5GaRWSmiKwVkX/12n+8iLwrIh+IyLWJ6lDVJ1X1e8B5wJmJynqN3+9JSDQLfWsOwFI4hNpDrAyttOkqckRaykZE7nGzLwmzMJGjo+soBx7AzAmNBs4WkdEiMkZEnu71il7VeKPzPUsRUiqOBJbiwW3SOYt3pGvZHBNjX0pOA6r6CvDvXru/DHygqstUtQP4PXCyqi5W1ZN6vdaK4R7gOZuSuripr6pneP1wq2gsCfHKGrGZa3NPSspGRC4VkcXAgY7Lc+S1HFjsgTxDgFVR263OvnhcARwNnC4ilySQe7KILBSRhevWrfNAzKwl6kw5U6fF4iWFPLTkpTWSjaFbr7PSlhqpWjaPAd8AnnLeI68vqep3PZBHYuzTeIVVdZqqfklVL1HVXycoN11VJ6jqhEGDBmUsZCRR5/LlJlGnVwonOlPnggULuO222+jq6vKmcoslCYU+tOSlNZKNodtCVtSFQErKRlVDqvoRcD7wFeC7wLnA5SLiRRjiVmBY1PZQ4BMP6vWUWIk6vWDatGncc889KWXqtFi8otCHllKxRtx0/F4O3Ra6oi4E0p2zeRI4GegCNke9MuV14AARGenkxzkLY0UVFJFEneBZos64mTrnzJnDc889h6pywQUXsGrVKs4//3xaW1u54IIL6OzszPzkFguF7xXo1hrJR8df6Iq6EEh3nc1QVT0+ebH4iMjvAD8wUERagVtUdYaIXA48D5QDM1X17UzOkw2ykKgzbqbOI444okdmzmHDhvXI3NmvX7/MT26xUByZUN0sJrUpywuTdJXNfBEZo6ppOwWo6tlx9j8LPJtuvbnC40SdcTN1Ro5FMnP2ztxpsXhJKUQGsCnLC5N0lc1E4HwRWQa0Yyb2VVXHeiZZHyNWps4Ikcyc8TJ3WiyFSj5i3uWr4y8FRZ1NbCDONPE611OsTJ2hUIgbbrihR2bO6MydFkshk8+Yd66G2wowhE4pk66yWYnxRNtXVW8XkeHAnoAN0JkmsTIV5iozp8WSDfpazDurvBKTrjfag4APiMy7tGHDxVgsligK3bvNklvStWwOU9VDRGQRgKp+5rgqWywWC2AnzS09SVfZdDpBMxVARAYBdvWhxWLpgZ00t0RIdxhtGvBnYLCI3AXMA6Z4JlUeUI0bFaeoKJXrsFgspUW6mTrniMgbQCPG7fkUVX3HU8lySHV1NRs2bKChoQGRWOHZigNVZcOGDVRXV+dbFIslJQo9Lbglc1JWNmJ646GquhRY6r1IuWfo0KG0trbiVUTofFJdXc3QoUPzLYalxPHS88qmBe8bpKxsVFVF5EngS96Lkx/69evHyJEj8y2GxdIn6Wsu0n2VdOdsFojIoZ5KYrFY+iTWRbpvkK432lHAxSKyAhPt2YarsVgsaWFdpPsGNlyNxWLJO9ZFuvSRdF1lRWQ34ABgu+uTqr7ikVxZQ0RCwPtRu+qBkMvtgcB6D8XpfS4vyscrE2u/m33R29lsi3jyZFI+0fF02qMv3Rux9pfSf8Xre6P3djG1RaIy0fv3UdX0Ux2rasov4CJgMfAZ8BKwFfhHOnXl+gVMT3cbWJhNWbwoH69MrP1u9vW6/qy1RTbaI9HxdNqjL90bya632NvD63sj0b1S6G2Rzr2RzitdB4EfAIcCK1T1KGA8UCx+w3/NcDubsnhRPl6ZWPvd7PtrgmNe43V7JDqeTnv0pXsj1v5S+q94fW/03i6mtkhUxrPrSGsYTUReV9VDRaQZEyetXUSaVXWcV4IVIiKyUFUn5FuOQsC2RU9se/TEtscObFsY0nUQaBWRXYEngRdF5DPgE6+EKmCm51uAAsK2RU9se/TEtscObFuQgYPA9gpEjsRMIv1NVTs8kcpisVgsJUValo2IVAPfx6SHVkwgznTnfywWi8VS4qQ7Z/M4JmHab51dZwO7qeq3PZTNYrFYLCVCusrmLVU9ONk+i8VisVgg/aGvRSJyeGRDRA4D/umNSBaLxWIpNdK1bN4BDgRWOruGA+9gsnWq9sEYaSJyEGb90UCgSVUfyrNIeUNETgFOBAYDD6jqC/mVKL+IyL7ADUC9qp6eb3lyjYj0Bx4EOoCAqs7Js0h5pa/eD+laNscDI4EjnddI4OvAScA3vBEtd4jITBFZKyL/6rX/eBF5V0Q+EJFrE9Whqu+o6iXAGUDR+tR71BZPqur3gPOAM7MobtbxqD2WqeqF2ZU0t6TYLqcCf3LuiW/mXNgckEp7lOL94IZ0lc0a4DTgF8DPMTfTGlVdoaorvBIuh8zCKNDtiEg58AAm6Oho4GwRGS0iY0Tk6V6vwc53vonxzGvKrfieMgsP2sLhRud7xcwsvGuPUmIWLtsFGAqscop151DGXDIL9+3RJ0l3UedsjDfar5zts4HfAEXpjaaqr4jIiF67vwx8oKrLAETk98DJqjoVY8HFqucp4CkReQZ4LIsiZw0v2sLJ5no38JyqvpllkbOKV/dGqZFKuwCtGIXTTIkukUixPZbkWLyCIN0f/kBVvVBVX3Jek4FRXgpWAAxhx9MYmD/MkHiFRcQvItNE5GHg2WwLl2NSagvgCuBo4HQRuSSbguWJVO+NBhH5NTBeRK7LtnB5JF67zAVOE5GHyH58vUIiZnv0ofuhB+laNotE5HBVXQAl640mMfbF9aZQ1QAQyJYweSbVtpgGTMueOHkn1fbYAJSi0u1NzHZR1c3A+bkWpgCI1x595X7oQbrK5jBgkoj08EYTkcWUjjdaKzAsansofSP+WyxsW/TEtkdsbLv0xLZHFOkqm+OTFyl6XgcOEJGRwMfAWcB38itS3rBt0RPbHrGx7dIT2x5RpDVn43icfQ7sAewTeRWrN5qI/A4IAgeKSKuIXKiqXcDlwPOYNUSPq+rb+ZQzF9i26Iltj9jYdumJbY/kpLuo8yLMAsaIh8nhQFBVv+apdBaLxWIpCdIdRotk6lygqkeJyH8At3knVvYYOHCgjhgxIq8yfLrpUz7+/GOzITCkbgh71u7Zo8zmjs0s3bAUFESEAxsOpH9l/53qclsuGyxZsoTu7m723Xdf+vePfc7mVc0AjBs2LmFd7777LgAHHnigJ+W85t0Nznkb4p/XTRnI3zXkA9dtkqdyFve88cYb61V1ULrfT1fZbFPVbSKCiFSp6lIRKYpfdcSIESxcuDCvMgRXBZn4vxMJa5iaihr+OOmP+Ib5epSZ+upUrv/H9QCUSRmTjprEdV/d2UvSbTnPryEYZOLEiYTDYVasWEFTUxM+n2+ncrtetSsAC+9L3OZ+vx+AQCDgSTmv8c9yznte/PO6KQPu26QUcNsm+SpncY+IZDRFku46m96ZOv9CH/aySBXfMB9j9xjLyF1H0jSpaSdFA+Af4adMzM9TWV6Jf4Q/Zl1uy3lNIBAgHA4D0NHRkfPOv5jpKu9iW/U2gquC+RbFYskZ6ToIfEtVN6rqrcBNwAzgFA/lKnnqq+oZXj88pqIBdwoplXJe4/f7KStzlFxl5XaLI9uEQiFWrlxJMFicHXVwVZDNdZtpr26ncXajVTiWPkO6w2jbUdWXvRDEsjOLLl7kqlx9VT31VfVJFY2XQws+n4+xY8cSCoWYM2dOzCE0rwkGg7S0tBAOh2lsbIw7dFfIBD4KmA8CHd0dBD4K5OwBwWLJJxkrm1Kgs7OT1tZWtm3blrNzPnSYyUDwzjvvZFzXLV+4pUdd1dXVDB06lH79+mVcdyLq6+upr69P2OF3re6iq7WLYDCYsWKINXQXr85CHbPfPsypUFmRu2HPUiHUHiK0LURwVTBjJV2o90ipkrKycXJTbFPVkone2traSl1dHSNGjMDEkCwuytab4awDBx6IqrJhwwZaW1sZOXJkXuUKBoNsnrsZuvHEEokM3YXD4ZwO3XmJb5iP/m396erXRdOFuRv2LAWCq4K0rGkhrGEaZzfmdNjYkjlJ52xEpExEviMiz4jIWmApsFpE3haR/xKRA7IvZnbZtm0bDQ0NRalowCiZAwcaZ0ARoaGhIadWWjwCgYAJKK/JnQgW7rGQBZ0LEs7FRIbuRo4cWZBDaKH2ECtDK5POw0x4YwKHLzjcdpQpEvgoQFgdy9YZgrQUD24cBF4C9gOuA/ZU1WGqOhj4KrAAuFtEzsmijDmhWBVNLOJdi9vO0Cv8fj+UA5LYiSBiAbUH22lsbEyocOrr6xk+fHjBKZrIU/fyjctzPvHvn+XfPiRUyuTL89LiDW6UzdGqeoeqtqg6jxWAqv5bVZ9Q1dOAP2RPRIsXZKMzbB7XTPO45rjHfT4f/U/tT5WvKqElkooFVKjYp+7sky/PS4s3JFU2qtoJICLfFpE65/NNIjJXRA6JLmMpXPLVGVbsVUH1odUJLRG3FlAhY5+6M8Ot1Z1syYClcEllnc1NqtomIhOBY4FHgYeyI1bh4/f7i6pTzEZnOG7cOMaNG5dxPW4toEKmrz11ux26c6NE8jkEackdqSibiPfZicBDqvoXoNJ7kSxz587lyiuv9LTOVDrDfMwBuLGAUiHX81Ngn7p741aJZMPqzsfvb0lMKsrmYyfl8ZnAsyJSleL3LS5ZtGgRhxxyiOf19pXOsNCflIs9CoJb3CoRr63uQv/9+yqpKIszMHkZjlPVjcDuwDXZEKoYyEaH8d5773HMMcdw3333ceutt3Lfffd5VncqFPJTYSAQSOpAUMiT9ZEoCMuXL0/qeeeWQv293CoRr4cg8/X79xWvwHRJRdlsBfoDZzvb/YCNXgtUDGSjw2hvb+eMM87g5z//OYMGDWLBggXcfvvtnq6XCZwXSLpaOpWnQjedXD6CThbyZL3XAUwL+Sk+FSXipdXt9e9vlYg3pKJsHsQkSYsomzbgAc8lKgKyEfH4xRdf5OCDD2bvvfdmwIAB7LnnnlRXV9PdndtADW6fCt10cvkKOpmNyXqvrAevA5gWshUH+Rm67WvOGsVCKsrmMFW9DNgGoKqf0UcdBLIR8bi5uZkxY8bw1ltvMXbsWNauXUtdXR1btmzh/PPPp7W1lQsuuIB169b12O7s9Nbr3O1ToZtOLvBRAIQeQSdzhdtOzo1XoZfWg9dREArZissnbn5/tw8QXg9T9lVLKRVl0yki5YACiMggIJz4K6VJNsKmDBgwgKVLl9LS0sLYsWO55ZZbuOyyyxg0aBDDhw/n6quvZtq0aTttex1s0+1ToZtOLpWO0Cs36mzgtfXgZRQE+xSfHm4fIAp5mLLYSEXZTAP+DAwWkbuAecAUL4URkZkislZE/hW1b3cReVFE3nfed4s6dp2IfCAi74rIcV7Kkgyvw6acc845vP/++9x+++089NBD7L777lxxxRVs2rSJZcuWUVFRQW1t7U7b2cDNU6GbTq5UOkLX1p6LObFs0Fe8DL3E7QNEoQ9TFhOulY2qzgF+AkwFVgOnqOofPZZnFnB8r33XAk2qegDQ5GwjIqOBs4AvON950LG8ipLdd9+dl19+mWHDhvHaa69x11130d3dzZVXXsmdd97JuHHj+Pvf/95jO1thXdx2mm46uVLoCEtFaVp2UOiZcEuRlFIMqOpSTNTnrKCqr4jIiF67Twb8zudHgQDwU2f/71W1HVguIh8AXwZyYudmo6Nvb2+nra2NhoYGACoqKpg5cyYA11xjvMyPPvroHtuW3OA2QZ2lOIg8QIS2hZhz6pykmXCTlbMkx7VlI4ZzRORmZ3u4iHw5e6JtZw9VXQ3gvA929g8BVkWVa3X27YSITBaRhSKycN26dVkVNhOqqqpYtmxZvsUoaPrq5KrFe9xa3V5b524dDsY/PJ59f7lvycwTper67KNwXJ9jxdHXWAVVdbqqTlDVCYMGDcqyWJZCIF/zJxZLIrx2TCimh69UhtEOU9VDRGQRGNdnEcmF6/MaEdlLVVeLyF7AWmd/KzAsqtxQ4JMcyGOxFBSloFSbm5vzLUJc3Lavm5TVsRwOYpV1W87LNNnZphhcn58CznU+nwv8JWr/WSJSJSIjgQOA13Igj8XBWg8Wi8G1JZKCY0JNRQ3lUh63XLG5ZRea6/PvMBP8B4pIq4hcCNwNHCMi7wPHONuo6tvA48AS4G/AZaqa2+X2FksRUWxpMYoJty7Sbj0bfcN8NE1q4o6j7ohbrtjcsl0No4nJM/wK8AbQiJkvOUVV3/FSGFU9O86hxjjl7wLu8lIGS34opuEAi6U3EYslrOGkLtJuPRt9w3yJU4GkcM5CwJWyUVUVkSdV9Utk0fW5kHn3XfN+4IH5laMUiQwHhDVM4+zGhE98VimlT6IU3hEilk8xpubOJ/lwkS42t+xUHAQWiMihqvp61qQpIiIeIHbOInPcToamopS8pvnu+8yH83JyuoLHKv2dycdarEUXL3JVrhD6q1TmbI7CKJwPRaRFRBaLSEu2BCs0uruhowM2bfKmviVLYPHi2PVlI1NnIeNl8E9LfLpWH8q2139IoowYodBoVq78TsIywVVB3lq4C8v/eib+O69LnIbig9GsfDpxfW5lA1j4XD8WPHJU0nKWwiMVy+aErElR4GzaBFu2mM/vvQejRnlfX3SYs2xl6ixU3A4HFNsYdSERDMLmx58FyjjqqG5eeqmc3mH9gkFobr4fEBoboamJncoAzH76fXTWi9BdScfLHcwe9yd8l8awRIPQPOV+UKHx2fj1BYOw+Y9/Ay1LeN5gEDbP/Qt0VyYsB0bJhZaOI3hgZmVKhUKwRFNRNqfF2BcSkTdUtdkjeQqStrYdn8Nhs53Jjxervtpak6nzsssuY8GCBTQ0NLBx40auuuoqby6iwHEzBJGNMerm5vtclevaugtdW2oJBouzY5o9ewUmwEY57e2dzJ7dis+3T48yZpomDFTQ3q4EAhL7Wj86ErorQSugW812DEx9CpTR0WG2Y9W3o5wkL9dVk7ScUXK/MsorjpILBqHlv35OuLNfQkUYKRsIgN+f+Ld3q7xyPZSVz+HnaFIZRpsAXIK5Y4cAkzExyx4RkZ94L1rhUFe343NZGSxpy8y/vXd9dXXxM3UuX77Ns6G7UiAfgT2DQdi8an/a1+9FYyNFOYTzafgPUNEB0gkVnWa7Fw0Ni4EOoJNweKuzvTOTTtkHqeiCsi6qqsqYdMo+Mcv5/VDWrxPKuqisNNvxylFuZEtarmJr0nKxlFesMtJdDVoetwyY33riRLj+ehL+9hHltXzuhUnvkWAQpk5Nfh+5KeemTOCjAOGVX4ZXr6X9o0PyNvycimXTAByiqpsAROQW4E/AERiX6J95L15hUFsLu+xi5m1GjoRfLXI3oe22vtpaePrpnpk6a2v3pLy8mtbWbj77bOehtmiWLOlZV6aMHw+hEMyZU5xP8anQ1bULXV2JLZZAAFCTBS7RE7XXeOkZtucRy6GzEVb6YXiAPY8Yt1OZDRuexqyV9lNW9iobNpwIjNmpnM8HB//kR4SWjmPO1ZPjtoXPBzWnfIOu1ok03XtLwnL9Tz3Zs3LblVe3UlnZL6ZS8vuhstLMwyZTXOIExkpmdYU7+/VQXvGsrokTzYhGTU3iIcPGxh3yxbPOkpUBaNhwEjz6A+iuJFzRQcOJH8a+2CyTimUzHPPYE6ET2EdVtwLtnkpVgIweDWPGmM7ci7Dj0fXBzpk6ly9fyy671PHii4/z6qvP8fnnygUXXMBDDz3Ec889h6o6mTu3smULtLeb+Z9EVtD48bDvvsmflFpaYPnyxE9ybutrvvu+HZ5cCVh4y3QWXPNY8glil+XcEAzC5s37096e2GLx+wFRQBN2TF4TCoVYuXIlwSQX62ax5qQjJsGwBfCVu6kauchsx6inrOw14G6qqt5MWGf9/ksYftJjSZVuxV6vU33oLzwrZ9Z8TyVRcHefD6qPOInyYbdz332LY9bp85nO+Y47Eg+hRZRSeXlipeTWioulvOKV6+jY4ZgUzzpLVgZgwztjtg97loVrzHYeSEXZPIbxRrvFsWr+CfxORPpjVvH3GbKR36R3ps5f/vIWvv3tyzjkkCNoaZnH3LkzOPPMMzn22GOZN28eM2aY7a6umu11ROZ/YuFWiQQCph5IfAO7ra9r6y5s2zA4qYJzM0zl9XCWubaeFkssfD7oP+wDqgauTtgxeUkwGKSlpYXly5fT2NiYVOEkwzfMR/9N/alqr+Kl816KnewuCxlo3dK1uottr29LeJ3BYJDNczfTHmxP2CbBYJBtr7xE96o7ueqqw+KWu+46P88/70/4e7pVSj4fjL3mR4z81kzPlFeycqnUVVNdQXk5VFVKzh6WepNK8rQ7gO8BG4EQcImq3q6qm1X1u1mSr2Dxeu6gd6bOPfbYnVtuuYLDDtuP1tZFLF3azHHHHcd+++3HokWLaG4223V1Zt4Hdsz/xMKtEvH7jXmf7AZ2U59b5RBrmCqTcm4x1+bOYqmo2UJ1w9qcDSsGAgHCTgN3dHQkHEr75JN9aGk5kenTY8+xRKjorqB6W3XiZHceZ6B1g1slEggEoBvQxG3itpxby9Hng+uuS/6Q4cbaS0V5JSvnZV25wPWcjROy5iCgXlVvj+SzUVUb/NIDIpk6Dz74YP7xj39sT6BWVwdVVRXcfPPN28tWVOzYrq018zltbaZsvDkbv98oo3A4sRKJ3JjJvG/c1Od2rmP7MJVCZYInL7fltpclsULy+aB//w/o6qqlqWnvgpqfMkNaZYTDYSorK+MOaU2fvpj3338IqOTiizuAxUyenJ9hklh0re6iq7WLYDAYV4HFUg6xyvr9figHuknYJm7KRSzHcDhMY2NjQksuGAwSCATw+/0JlbBbLzOfz12H76acl3Vlm1QcBB7E+EV+Dbgdk8/mCeDQLMhV8GTDfbF3ps5QKMQNN9zAueeey+DBg3fajlBbm9wxwOeDefPcuXC6vcmT1edWOfh8MO76y11NOM//Z5mra3BLRcUWKiq24PPtnXllHhIZ0gqFQsyZMyduJ/fEExswz4AVgPLEExuYPDl2neOax2VJ2thELBa6Sdihu1UiPp+P/qf2p6u1i6Z74ysHN+ViWY6xygWDQSZOnEg4HKampibnw4ulRDHks+kz9M7UWV9fz/333x93O1W8frpJVl9krqNrSy1NTyW2HOr3X0L9/kvw+eL0lC7PGSEUMq9k62LGjbvK+RRIWN+4a5OX8zquWH19PfX19Qk7t9NOa+CFFzoww4GdnHZaQ9yyuY535tZicatEACr2qqBir4qkHX6ycm4tR7dKyZKcVJRNvvLZWIqYCbcZ5eHzBXJ2zojzQjhM0pXmXhIKhQiFQgmHjFLBjXKYPHkM9957LuvXf5G77z4+4yE0twrJTbIztxYLuFciXuHWcvT7/dTU1NDR0ZH0GiyJSUXZRPLZ7OHkszkduDErUlksGRDLeSHbfVgqcwBe8957j+bkPNG4mYtJxWLxmnHjxiUt48Zy9Pl8NDU1uZqzcTu347ZcqeFa2ajqHBGJ5LOBLOSzsVi8wK0zBLgLu++GbAy3uA2TkmvczsVA7i2WVHBrxfl8vqTyB4NBGhsbt1tA8drEbblSJKmyEZEfxTl0goicoKo/91gmiyUjfD4YOza3URDczgGAu3hsblea5wO3czHZwI3Fkg8CgQBbt24FEreJ23KliJt1NnXOawJwKTtio10CjM6eaBZL+ixf3sy//92csw7a6wWRbtdFeY1/ln977pO4ZSJzMZJ8LqavEJnbKS8vT+qW7aZcKZLUslHV2wBE5AVMbLQ2Z/tW4I9Zlc4FInI88EvM7f8/qnp3nkWyFABuvcy85P33/0BXVy2QuRt1ZHFtsthd+SCVuZhCtUS8xu3cTipzQKVGKg4CvWOjdQAjPJUmRRzvuAeAY4BW4HUReUpV+1T4HEv+CQZh69ZRnnnAuV1cmy8mnDABoE91lslwM7eTSrlSIxVl8xvgNRH5M8b9+VtA7t1gevJl4ANVXQYgIr8HTqaPxWqz5J9seMAVwqpvi8UrUvFGu0tEngO+6uw6X1XdJcDOHkOAVVHbrcBhOTmzm3golpIh2bqSVDzgLJa+iBtvNFFVBVDVN4E3E5XJMRJj305yiMhkTLI3hg8fnm2ZMmbu3LkEAgGmTZuWb1FyRq6zF3pNKh5wbnLoWCylhhtvtJdE5AoR6dFLi0iliHxNRB4Fzs2OeElpBYZFbQ8FPuldSFWnq+oEVZ0waNCgnAmXLosWLeKQQw7JtxiWFFm0CJYtS6xA3ObQsViCwSBTp07NOL1EoeBmGO144AJM7pqRmBQDNRhF9QLwC1VtzpaASXgdOMCR62PgLOA7OTmz2+BbKfDee+9x2WWXsWDBAhoaGti4cSNXXXWVJ3VbCoNYOXSK0bopdku00HEbALSYohG4cX3ehon4/KCI9AMGAltVdWOWZUuKqnaJyOXA8xjX55mq+nbWT5yF4Fvt7e2cccYZ/OY3v+Hkk09m/vz5jB49mksuuYTq6mqPBLfkmx05dJKnSehL9BUXabe4iUhRbBGpU8nUiap2qurqQlA0EVT1WVUdpar7qepdOTlpFlbcvfjiixx88MHsvffeDBgwgD333JPq6mq6u7szrttSOERy6FRV5S7rp6X4cLP4M5UEe4UwJJeSsrE4RFyPwDPXo+bmZsaMGcNbb73F2LFjWbt2LXV1dWzZsoXzzz+f1tZWLrjgAh555BGee+45VJULLrhge+gLS/FQUbGF6urcZf20FB+RxZ933HFHwjxAbqIRRCyg66+/3pMU4+mSyjobS4QsBN8aMGAALS0tVFRUMHbsWG655RYuu+wyBg0axPDhw7n66quZMWMGa9asYebMmXz88ceceeaZ1NTUeHBBFoul0Ei2+NNtNIJAIIBJtJzfeGwpKRsnWVq5qtrH6fp68/LoRzvnnHP41re+xdy5c9ltt90466yzuOKKK9i0aRPLli2joqKC2tpaamtrWbRoEaFQiIsuusiTc1ssxYJ1TOiJm2gEfr+fysrKvOfkcT2MJiI/AFYDH4jIO87EvMUjdt99d15++WWGDRvGa6+9xl133UV3dzdXXnkld955J+PGjds+JltRUcHNN9+cX4EtFktR4GZILhe4WdR5H2Yh5w+Ag1R1rZOl8zYRuUNVb8qyjIVJFiIHtLe309bWRkODSe1bUVHBzJkzAbjmmmsIhUJcfvnlnHvuuQwePNjz81tyQz6ChFr6NoUQj83NMNrLwHiMy/N8EfkcaAEWA5eIyH8XkndaMVNVVcWyZcviHq+vr+f+++/PoUQWi8XiDW7W2fwZ+LOIHA78EDOUdjAwFtgdCIhIrarun1VJLRaLxVK0pOIgcBnwONCMsWoOAharqt9xHLBYLJYeuJ3QtxP/pU8qUZ/fF5HDMLljxmGG0n7iHOtI8FWLxWLJKVZ5FR4puT47SuUZ52WxWNIg0Upvi6VUsREELH2artVdbHt9W8lE1rVYChUbQcCSVQp5OCMYDLJ57mbohsbGxoIPZGixFDOpLOoUETlHRG52toeLyJezJ5rFkl0CgQB0A5o8kOG45nGMax6XI8ksltIjlWG0BwEfcLaz3QY84LlERYLfb1P/Fjt+v98kphDyGsbDYukLpKJsDlPVy4BtAKr6GWBdnrPA3LlzufLKK/MtRsnj8/nof2p/qnxVdgjNYskyqczZdIpIOU7mJydkTTgrUvVxbFro3FGxVwUVe1VYRWOxZJlUlM004M/AYBG5CzgduDErUhUBWcgKbdNCW3oSGdazrtKWEsD1MJqqzsEs4pyKCVlziqr+MVuCFTKRrNDLl5us0F54zUbSQv/85z9n0KBBLFiwgNtvv51t27ZlXrnFYrHkmVQXdS4FlmZJlqIhVlboTK0bmxbaYrGUMqm4Pn9bROqczzeKyFwR8WRiwan7bREJi8iEXseuE5EPRORdETkuav+XRGSxc2yaRFLR5YAsZIWOmxb68ccf75EGetWqVT3SRHd2dmZ+ckthEgrBypXemM4WS55JxRvtJlVtE5GJwHHAo8BDHsnxL+BU4JXonSIyGjgL+AJwPPCg46SAc+7JwAHO63iPZElKJCv0yJHQ1OTNnM2AAQNYunQpLS0tPdJCH3HEEcybN48ZM2Zw5plnMmzYsO1poqdNm0a/fv0yP7ml8MjGWK3FkkdSUTaR8ZwTgYdU9S945Pqsqu+o6rsxDp0M/F5V21V1OfAB8GUR2QsYoKpBVVVgNnCKF7K4pb4ehg/3zjngnHPO4f333+f222/noYceYvfdd+eKK65gv/32Y9GiRTQ3N3PcccftlCbaUqLEGquNh130ZSkCUlE2H4vIw8AZwLMiUpXi99NhCLAqarvV2TfE+dx7f0xEZLKILBSRhevWrcuKoJkSKy10ZGQwkga6q6srZppoSwmSjbFaiyWPpOIgcAZmqOpeVd3oWBfXuP2yiPwd2DPGoRscKynm12Ls0wT7Y6Kq04HpABMmTIhbLhWy0c/3TgsdCoW44YYbeqSBjk4TbSlhImO1oRDMmZPYhHbrh29dqS15JJV8NltE5CXgABE5wtnt2i9XVY9OVTiMxTIsanso8Imzf2iM/UVN77TQNg10H6e+3rwSKZDI3E44bOZ2Ek0iZmNxmMXiklS80S7CTOA/D9zmvN+aHbG28xRwlohUichIjCPAa6q6GmgTkcMdL7RJQDzryGIpTgKB5FaI27kdlw4H029ZyGPXLPDOISEYhKlTk9eXj3Ju67J4QirDaD8ADgUWqOpRIvIfGKWTMSLyLeBXwCDgGRFpVtXjVPVtEXkcWAJ0AZepasRR4VJgFlADPOe8LJa+RWRuJxxOPLfjZnFYMMj+qzYjSnIrKRg0dfj9ics0NprzVVbGry8YhIkTjXw1Nd6US3Zet7K5vdZslCsxUlE221R1m4ggIlWqulREDvRCCFX9MyYUTqxjdwF3xdi/EPiiF+e3WIoWt3M7bpRSIICoMyGaaLWy204/EICtW83nRPUFAhBZJudVuY4O6O6OX86tbKkoTC/LRcp6pZQKQMGlomxaRWRX4EngRRH5jBKYJ4mgquRwXWhWMd7glpzjdgLe64l6N3M7Ph/Mm5e4w/H7aa8so6IrTEWmVpJTHzU1OzrWePX5/eZ4Lsu5lc2N4kq1nNdKzisLM8uk4iDwLefjrY6jQD0lMnRVXV3Nhg0baGhoKHqFo6ps2LCB6urqfItiyRVulZbPl1Qh/eiasYxbGmLy1UmsJDcdtc9nOrZknWE+yrmtKxuK0Csl51aJuFWEWca1snHW1ZwGjIj63jjgds+lyjFDhw6ltbWVRGtwVv57JQDDdx+eK7HSprq6mqFDhyYvaHHNfc3N3lVWwF5hS/avZ8n+9UxOZiW56agjZd1cYz7KuS2TD4XpRnm5VSJuFWGWSWUY7S9ACHgDaM+OOPmhX79+jBw5MmEZ31XmR9x438YcSGQpWVJxVS5k3Hb6pUC+FGEypeRWiaTycJBFUlE2Q1U1Z/HHLJaiw43Fko2Q4ZbSxMWwp+cWZhZJRdnMF5Exqro4a9L0ISL57hOFnHFTxlIguLVY3LoqWyxuKAAl4pakizqdMP4twETgTSfUf0vU/j5B1+outr2+jaBHC8BGh0J8J0n4+H0++YQTW1pYPH26J+f0+/3bFVjB4fECu/uam72dZ0mG28WV2QgZbrEUAW4sm5OyLkWBEwwGGfPKZvzV8NNvH8E9f3wls5z1wSA/f+stKlTpPuooyl96aadOZ/H06Tz0/vtUAh0XX8xiYMzkyXHry/d4bFzGj0++BsTtuo0U2KWri9quLs8m4ZPWl4rFsmhRxvJYLMVGUstGVVeo6gpgDcYb7RfAzzH5Z9ZkV7zC4NVZP6NpPdyxAv62potXZ/0so/pWzJ5NhSoVQLi9nRWzZ+9UZsMTT1CJeRro52zHJOL+eNNNSfOeuLaU3FoZ48fDvvvGL+c2J0sq4fSTndM57/6bN7NXe3vSNukq72Jb9TaCqzKsz1osFktCUkkRMBuTxOxXwP3AQcBvsiFUoTFq4ydUdkOFQr9us50JLwMdQKfzejlGmYbTTutRpuG002JXFlkkFu3+GIOIpfTDzz5jv4svjq9w3CovN4rErRKJrD0oL09sFaSgvIReK+FjVbcqSNcum+msaqdxdmN8heOyPs+THFksJUQqyuZAVb1QVV9yXpOBUdkSrJDY7+wL6SiHToHOcrOdCQdMmsQxItwMfL2ykgMmTdqpzJjJk7n0gAO4b7fd+PDhh+MPofn9bCsrowsSdtSuLaVAgC5HeWl7e/yO1Y0icZuTJeJVc8cdia2CFJRXNybbX3dFRdzzBj4K0FkG4TLo6O4g8FH8+hQnh0UOJ/X9s/z4Z+XmXBFC7SFWhlYmtvQsljRIxRttkYgcrqoLAETkMOCf2RGrsBhzymS+dvoP+MonnZzxwwcZc0qcjt8lPp+PLQcfzB9CIebMmRN3/ufR995zUxkX77cfX1y/nuPvvpsxcepqOO00Ol54ASWxpbS4oYH9MB1rZzjMhw0NjIlV0M1KaDchUqLLJrMIXK6+DgLXAl8B5qsyFYhVs3+En35h0x6VlZX4R8SuD5+PD/r3p7ari72TKcMiJrgqSMuaFsIapnF2I02TmvANs1aaxRtSUTaHAZNEZKWzPRx4R0QWA6qqYz2XroB4c48q3tyjijsyVDQR6uvrqa+vz8zRAOO88NsPPyQcDnPLVVfRNGZMzDrHTJ7Muffeu0MpxbGUnt6wgacAP/BqWRknbtgQW9mksmLaq2Ell+cMBAK8gsmHUd7dTSAQiNkmvmE+fve4AMr+99zHmAQd65aKCrZUVLB3DofIQu0hQttCBFcFc9LpBz4KEFZjOUYsPatsLF6RirLp0ws6x40b52l9Xq2dCQQChJ2hpY6OjrgdK8CKvfdmxd57c028ITmMe/SNZWUsCIepqarivxINGeXDx9/FOf1+P2VlZYTDYWOxJJgD+ub7SpmCfOcqaBpTMPMt+bAy/CP8lEkZYQ1TWZ7A0rNY0sD1nE3EKy3eK5tCWuIT6ViBxB0rRjElU3I+n4+xY8cycuRImpqaMra88oHrawgEjKKBpF5wb3V1MWObd+uskhHLysg2vmE+xu4xlpG7jrRDaBbPScWy6dMEzgvkW4SYRDrWUJL5n1RYVALrQFwNU/r9hAVEoSzRHFAwyAVbtxIOh5nR2JgTJZwvK6O+qp76qnqraCyeI30t94mIhID3o3bVYwKMutkeCKz3UJze5/KifLwysfa72Re9nc22iCdPJuUTHa8HQnXQfwDUfQ5tbebhK1Z71ABDnG3F5HH6lGzfG5V0UUUd7bTRweak5b29N2LtL6X/StJ7I419ufqv5Kvf2EdVB6Vw3p6oap96AdPT3QYWZlMWL8rHKxNrv5t9va4/a22RjfZIdDyd9uhL90ay6y329vD63kh0rxR6W6Rzb6TzSmWdTanw1wy3symLF+XjlYm1382+vyY45jVet0ei4+m0R1+6N2LtL6X/itf3Ru/tYmqLRGU8u44+N4yWCSKyUFUn5FuOQsC2RU9se/TEtscObFsY+qJlkwnehF8uDWxb9MS2R09se+zAtgXWsrFYLBZLDrCWjcVisViyjlU2FovFYsk6VtlYLBaLJetYZeMRInKQiPxaRP4kIpfmW558IiKniMgjIvIXETk23/LkGxHZV0RmiMif8i1LPhCR/iLyqHNPfDff8uSbvno/WGUDiMhMEVkrIv/qtf94EXlXRD4QkWsT1aGq76jqJcAZQNG6OXrUFk+q6veA84Azsyhu1vGoPZapamZJkAqMFNvlVOBPzj3xzZwLmwNSaY9SvB/cYJWNYRa9olqLSDnwAHACMBo4W0RGi8gYEXm612uw851vAvOAptyK7ymz8KAtHG50vlfMzMK79iglZuGyXYChwCqnWHcOZcwls3DfHn0SG4gTUNVXRGREr91fBj5Q1WUAIvJ74GRVnQqcFKeep4CnROQZ4LEsipw1vGgLERHgbuA5VX0zyyJnFa/ujVIjlXYBWjEKp5kSfcBNsT2W5Fi8gqAkf3iPGMKOpzEwf5ghccoiIn4RmSYiDwPPZlu4HJNSWwBXAEcDp4vIJdkULE+kem80iMivgfEicl22hcsj8dplLnCaiDxE9kMeFRIx26MP3Q89sJZNfCTGvrgrYFU1AASyJUyeSbUtpgHTsidO3km1PTYApah0exOzXVR1M3B+roUpAOK1R1+5H3pgLZv4tALDoraHYkLL90VsW/TEtkdsbLv0xLZHFFbZxOd14AARGSkilcBZwFN5lilf2LboiW2P2Nh26YltjyissgFE5HdAEDhQRFpF5EJV7QIuB54H3gEeV9W38ylnLrBt0RPbHrGx7dIT2x7JsYE4LRaLxZJ1rGVjsVgslqxjlY3FYrFYso5VNhaLxWLJOlbZWCwWiyXrWGVjsVgslqxjlY3FYrFYso5VNhZLDhCRK0XkHRGZk29ZLJZ8YNfZWCw5QESWAieo6vKofRXOwj+LpeSxlo3FkmWcCL/7YtJPhERkuoi8AMwWkUEi8oSIvO68vuJ8p0FEXhCRRSLysIisEJGBeb0QiyUDrGVjseQAEfkIk8H1cuAbwERV3SoijwEPquo8ERkOPK+qB4nINGC9qt4uIicCTwODVHV9vq7BYskEm2LAYsk9T6nqVufz0cBok28OgAEiUgccgUmnjKo+IyKf5V5Mi8U7rLKxWHLP5qjPZYAvSvkA4CgfO+xgKRnsnI3Fkl9ewAytASAi45yPrwDfdfadAOyWc8ksFg+xysZiyS9XAhNEpEVElrAjg+NtwBEi8iZwLLAyXwJaLF5gHQQsliIg4mBgHQQsxYq1bCwWi8WSdaxlY7FYLJasYy0bi8VisWQdq2wsFovFknWssrFYLBZL1rHKxmKxWCxZxyobi8VisWQdq2wsFovFknX+P7q4WmEOU8TMAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "res = ImpME\n",
    "rho = 1e12 * np.abs(res.impedance)**2 / freq[:, None, None]\n",
    "rho_err = 1e12 * np.abs(res.error)**2 / freq[:, None, None]\n",
    "phi = np.angle(res.impedance, deg=True)\n",
    "rad_err = np.arcsin(res.error/abs(res.impedance))\n",
    "rad_err[np.isnan(rad_err)] = np.pi\n",
    "phi_err = np.rad2deg(rad_err)\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = plt.subplot(2, 1, 1)\n",
    "ax.set_xscale(\"log\", nonpositive='clip')\n",
    "ax.set_yscale(\"log\", nonpositive='clip')\n",
    "ax.errorbar(freq, rho[:,0,0], yerr=rho_err[:,0,0], fmt='k.', label=r'$\\rho_{xx}$')\n",
    "ax.errorbar(freq, rho[:,1,1], yerr=rho_err[:,1,1], fmt='g.', label=r'$\\rho_{yy}$')\n",
    "ax.errorbar(freq, rho[:,0,1], yerr=rho_err[:,0,1], fmt='r.', label=r'$\\rho_{xy}$')\n",
    "ax.errorbar(freq, rho[:,1,0], yerr=rho_err[:,1,0], fmt='b.', label=r'$\\rho_{yx}$')\n",
    "plt.xlabel('freq')\n",
    "plt.ylabel(r'apparent resistivity  $\\rho$ ($\\Omega.m$)');\n",
    "plt.legend()\n",
    "\n",
    "plt.title('Site 002 results in 2-stage RR ME\\n  configuration with sites 99 and 100 as remote references')\n",
    "ax = plt.subplot(2, 1, 2)\n",
    "ax.set_xscale(\"log\", nonpositive='clip')\n",
    "ax.errorbar(freq, phi[:,0,0], yerr=phi_err[:,0,0], fmt='k.', label=r'$\\phi_{xx}$')\n",
    "ax.errorbar(freq, phi[:,1,1], yerr=phi_err[:,1,1], fmt='g.', label=r'$\\phi_{yy}$')\n",
    "ax.errorbar(freq, phi[:,0,1], yerr=phi_err[:,0,1], fmt='r.', label=r'$\\phi_{xy}$')\n",
    "ax.errorbar(freq, phi[:,1,0], yerr=phi_err[:,1,0], fmt='b.', label=r'$\\phi_{yx}$')\n",
    "plt.xlabel('freq')\n",
    "plt.ylabel(r'phase $\\phi$ (degrees)');\n",
    "plt.legend()\n",
    "plt.ylim(-180, 180)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now computing 2-stage Bounded Influence Transfer Function for site004 with sites 100 and 99 as remote references for a rejection percentage of 1% and 3 bounded influence steps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SignalSet: 9 channels, 1 run\n",
      "tags: {'site004_Ex': (0,), 'site004_Ey': (1,), \n",
      "       'site004_Hx': (2,), 'site004_Hy': (3,), \n",
      "       'site004_Hz': (4,), 'site099_Hx': (5,), \n",
      "       'site099_Hy': (6,), 'site100_Hx': (7,), \n",
      "       'site100_Hy': (8,), 'B': (2, 3), \n",
      "       'E': (0, 1), 'Bremote': (8, 6, 5, 7)}\n",
      "----------  -------------------  -------------------\n",
      "  sampling                start                 stop\n",
      "       128  2016-04-28 19:00:00  2016-04-29 07:00:00\n",
      "----------  -------------------  -------------------\n",
      "starting frequency 0.001\n",
      "starting frequency 0.00139742\n",
      "failed to converge (maxit=100). while processing step 2 (weighting=1).\n",
      "starting frequency 0.00195279\n",
      "starting frequency 0.00272887\n",
      "starting frequency 0.00381338\n",
      "starting frequency 0.00532889\n",
      "starting frequency 0.00744671\n",
      "starting frequency 0.0104062\n",
      "starting frequency 0.0145418\n",
      "starting frequency 0.0203211\n",
      "starting frequency 0.0283971\n",
      "starting frequency 0.0396827\n",
      "starting frequency 0.0554535\n",
      "starting frequency 0.0774919\n",
      "starting frequency 0.108289\n",
      "starting frequency 0.151325\n",
      "starting frequency 0.211465\n",
      "starting frequency 0.295506\n",
      "failed to converge (maxit=100). while processing step 4 (weighting=1).\n",
      "starting frequency 0.412946\n",
      "starting frequency 0.57706\n",
      "starting frequency 0.806396\n",
      "starting frequency 1.12688\n",
      "starting frequency 1.57472\n",
      "starting frequency 2.20055\n",
      "starting frequency 3.07509\n",
      "starting frequency 4.2972\n",
      "starting frequency 6.005\n",
      "starting frequency 8.39151\n",
      "starting frequency 11.7265\n",
      "starting frequency 16.3868\n",
      "starting frequency 22.8993\n",
      "starting frequency 32\n",
      "(32, 2, 2)\n"
     ]
    }
   ],
   "source": [
    "from razorback.weights import bi_weights\n",
    "from razorback.prefilters import cod_filter\n",
    "\n",
    "sig = prepare_signalset(inv, 'site004', ['site100', 'site099'])\n",
    "print(sig)\n",
    "ImpBI = rb.utils.impedance(\n",
    "    sig, freq,\n",
    "     weights= bi_weights(0.01, 3),  # bounded influence with reject probability of 1% and 3 steps\n",
    "    remote='Bremote', # including the remotes references in the computation,\n",
    "    prefilter=cod_filter(0.0), # prefilter: cod_filter(0.0)\n",
    "    fourier_opts=dict( Nper= 8,  overlap= 0.71) # fourier options with 8 periods by window, and 71% of overlap\n",
    ")\n",
    "print(ImpBI.impedance.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-19-5eea7cae3846>:5: RuntimeWarning: invalid value encountered in arcsin\n",
      "  rad_err = np.arcsin(res.error/abs(res.impedance))\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(-180.0, 180.0)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZsAAAEpCAYAAABFmo+GAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABkOUlEQVR4nO2deZwcVbX4v2dmMgvJpIFJWJMhYcsDSUgAgdYIjYMsig8UIaA+diPIIu+HqIjskKDPpxg2AxIhGlQ2lYcg4EADMY0sZhg0JAGSkExYAiNpZkJm7fP741ZPejq9VHdXr3O/n09/Zqrq1q1Tt6ruucu554iqYrFYLBZLPqkqtgAWi8ViqXyssrFYLBZL3rHKxmKxWCx5xyobi8ViseQdq2wsFovFknessrFYLBZL3rHKJkNE5Gsi8kSx5ShVRGSSiKiI1HiQ12MicroXclksluJilU0CRGSmiCwRkbCI/FtE/iYinwRQ1UWqelRMWhWRPXO41iQReVpEPhaR5SJyZNzxr4rIWyKySUT+KCLbxxz7iYi8LiJdzrmnZStHvhCRoIick825qnqsqt6TxTXrROQup9y6RGSpiBybjQwicrWI/Cabc3NFRM4QkUER6RaRj0TkFRE5LuZ4VLF3O781IvL9NHmq8y51i8gHIvJbEdk25riIyKXOe7VZRNaKyI0iUheTZqtnKiIBEemIu86rIlIVs+96Ebk7iezvicgjIvK5uHzXOHJ0x/xuSVA+3SKyWkR+JSJ7p7j/eDmDItIjIhNj9h0pImvSyLBLsoZV/DUSlZvzXvXH5bnRbfk527VOPq87z3SNiCwQkUlx9xZ7jf9LVjb5xiqbOERkLPAIcDOwPbArcA3Qm6dL/hZYCjQBlwMPiMh4R5ZPAPOB/wJ2BD4Gbos5dxPwRcAHnA78XEQ+5eai8R9IhVEDrAMOx5TNFcB90Y+wzAip6hhgW8yz/12scnDY1knzFeCK+Ao7Afs76XcHtgOujjk2D5gNnAY0AscCnwXuy0L2XYBT0qSJyr4/8CTwBxE5Iy7NF1V1TMzvgphj0fLxAUcCm4GXRWS/DOTchHlHUhEvw9sZ5J+M38fluW3c8XTl9wDwn8BXMfe/P/Ay0BKT5oK4a3zRA7mzQ1XtL+YHHARsTHH8DGCx8/+zgGJe1m5glrP/OKAN2AgsAaYlyWtvjBJrjNn3HHCu8/8c4N6YY3sAfbHp4/J7GLgkybEA0AF8D3gX+DWmsfF94E2gE1OhbO+krwd+4+zfCLwI7OgcWwMcGZP31cBvnP8nOWVSA9wADAI9TvncAgjwM2ADEAbagf2SyBwEzoktd+AnwIfAauDYDJ5rO3BiiuPfA9YDXcAKzAd7jFPe/Y78rzhpzwRec9KuAr4Zl9d3gXeAt4FznPLY0zlW59zDWuA94BdAQ7p3zdnexsnrk/FlHZPmBeDSFPc5JIuz/S3gCef/vZzndXDcORMx7+ln459L/PsVd53vAa9H5QOuB+5OJruz/ztOuVQletdSlU/M/keAB1J9B3Hv2FXOs4w+oyOBNTFpEsqQ4h6GXSPJ+3w1zjeT4jmlKr+oYp2YIo+tnlMxf7ZnszUrgUERuUdEjhWR7ZIlVNXDnH/3V9Nq+L2IHAAsAL6J6a3MBx6OHYaI4RPAKlXtitn3irM/evyVmOu9ian8thomEJEG4JPAv1Lc206Y3tpumNbrRcAJmB7ALphK/FYn7emY1tJE5z7OxbzcrlHVyzHKM9q6ugA4CjjMuYdtgVkYheaGQzCKYBzwY+AuEZF0J4nIjs71EpaNiEwBLsBU4o3A0ZjK5i8YhR9tge7vnLIB06AYi1E8P3OeOyJyDPD/MJXBnpiyjeVHjizTneO7Ale6uIdq51r9wFtJ0hwK7Ae8kS4/J/12mOf/vLOrBVNJvhCbTlXXOWnS9ZjieQj4CKMUMjlnB2BKhteKz+MzGaRfD9zJ8B5eKZCq/I4EXnCeTVlglU0cqvoRMBPTsrgTeF9EHnYqLDd8A5ivqn9X1UE1cw69wKEJ0o7BtO5jCWOGL9wcj+UXGMX0eArZIsBVqtqrqpsxCvFyVe1Q1V7Mx/YVZ4itH6Nk9nTu42WnbHKl35H/PwBR1ddU9R2X576lqneq6iBwD7AzZngxKSIyClgE3KOqy5MkG8T0OPYVkVGqusZR7AlR1T+r6ptqeAZ4gi2V28nAr1T1X6r6MWYINiqLYN6P/1bVfzuNjDmkHio51BnL78H0iL6uqhvi0nwgIpuBEGao7Y8p8gP4h5PnB0AzpkEERoknexbvOMczQTHDU1cmaWwlIjo8tX3Mvj+KyMaY3zdc5LF9mjTxzAW+6AxdJyJWhj9mmHcyTo67r6fjjqcqvyaSP6tY5sVd4zovBM8Gq2wS4FSAZ6jqBExLcRfgJpen7wZcEvuAMb2DXRKk7ca0jmMZi+nSuzkOgIj8jyPnyer0n5Pwvqr2xMn6hxg5X8NUvDtihtkex8wRvC0iP3Yq7pxQ1acww2m3Au+JyB3OPJkb3o3J52Pn3zHJEjuTq7/G9AYviNn/WMyE6ddU9Q3gYoyy3SAivxORRM8rev6xIvK8GOORjcDn2VIR74KZL4oS+/94zFDYyzFl/hdnfzKeVzOWvx1mmDRRi30cphy+gxnCSfecDnDyrAduB54TkXqM8tk5yTk7O8cBBhJcYxSmITEMVX0UM2Q4O41MUXZ1/v47Zt8JqrptzO9OF3n8O02aeDnfx7yX1yZJEivDCWmyS1Q+sHUZ3Rd3X0ckkCtZ+XWS/FnFclHcNdLNTeUNq2zS4LSG78ZU5m5YB9wQ94C3UdXfJkj7L2B3EYntqezPluGefznbAIjI7pgW+MqYfddgJnGPctHziFdE6zDzHrGy1qvqelXtV9VrVHVf4FOYYaOotdsmTKUZZacMromqzlPVAzHDhHsDl6aRO2OcXsRdGMV5oqoOfeRqrNyiE6aLnH33qupMjAJWzHDXVvI7LcwHMb2MHZ1K+1HMXBSY1uaEmFMmxvz/AWYo8hMx5e1TM8GdElXtxsyv/JeIzEhwfFBV/xfTA/pWuvycc/qBXwKTMe/3U8BEETk4Np1jqXUo0OrsWouZr4hlMkmG94AfYoxftklyPJYvYYYpV7hImyqP57I473+AI4ADc7g2mPIZJyJDz9V5H3cjeRmlIlH5/RU4WEQmJD6l9LDKJg4R+Q8RuST6EJ0P7VS2jGvH8x7GqifKncC5InKIGEaLyBfiFAoAqroSY0hwlYjUi8iXgGmYygzM8M8XReQzIjIa0+p6KDrHIyKXYSxRPqeqbuc9YvkFcIOI7ObkN15Ejnf+P0JEpjpzBR9hWmSDznltwCkiMkpEDsJYQSVjWPmIyCedshmFUVo9Mfl6ye3APhgropRzTSIyRUQ+6yiSHoxCiMr0HjBJtpig1mIU/vvAgBiT6qNisrsPOFNE9hGRbYiZj1HVCOb9+JmI7OBce1cROdrNDTnP+JeknuO5Efiu01NJScw80GbM3OFKzDuxSEQOFZFqZ1jpQeCvqvpX59TfO/d4sPOO7w38N/C7JHIHgVcx84DJZNlRRC7ATNZf5pSVaxxZJ4vIzZje3TVpTkkk50bgfzEGHplQ53y/9U65dwB/B34kImOc9+pSTI8nWT2SSq4gceXnPIuo9d6BIlIjIo0icq6InJXpNQpCvMXASP9huuD3YSYNNzl/5wNjneNnMNxC6FxMa3YjZhgLjBXTi86+d4D7SW5BNgljNbIZ05o7Mu74VzEtpU3An3CsxXSLxUovZrgt+vtBkusEiLOQwTQ2/p9z3S6MVdoc59ipzv5NmAp3HlusYnbHfEzdwJ+dY1tZoznbfkxP7EMnXQvGMqwb09JfBIxJInOQOGu0uOPDLKti9kd7J1EruOjva0muMw1jxdWFGX55BNjFOdaEsYL7EPiHs+98p0w2YobpfgdcH5PfZZghv7eB8xxZJjrH6jHzNKswSvw1zFBHIrkS3fME55lPiy9r57hgesQXJskz1nryI8x7enTcO/E9jJHBZkzv98dAfVw+ZznX+chJ+30cC7JEzwZj3KFsbY3W7cizAdNDPCbuOmscOWKf4x9iymcwJo+3MHN5+6T4vgNsbY12Tsz2GEeWNXEypLJGi/8dienR3u+8Bx9ghqT3jTn3arZYOcb+dnBTfs6+WoxSfSPm/n8JNMfcW/w38HIh6tFEP3GEygqntd2jZsLWYrHEISL7AP8E6lR1oNjyWCzFIqNhNBGpErOi/c8isgFYDrwjIv8Skf8Rkb3yI6bFUj6IyJfErO7eDjP3839W0VhGOpnO2TyNWVh4GbCTqk5U1R0wFjLPAzeKyNc9ltFiKTe+iZnTeRMzzHNeccWxWIpPRsNoYtYgbGXemGkai8VisYwssp6zsUrFYrFYLG7JyvRZRH6JWZC3TkT+LiJ3isiFHstWFjgmm8+K8S78vyLyA6d8SgIR+YWIFG0hV7rri0delYt9n8VGjAfkxcWWw5IdYjw6fyAi76ZPXZ5ku87mM5gFbROBLwN/AEZ7JlV5MRtj2jhWVS9R1TmqmpVL/VxJVOGo6rmqWjQXFbHXlySu18vpOs76mafEhJ94w1kbFXv8HGd/t4j8RVJ4IigkYjw1rBCRiGztVRkR+W8Rede5rwUyPKTA9iLyBzFu7N8Ska8WVPgC41UDKIPrTQQuwZhGp1ogXdZkq2yex7jPQM1q80dV9UbvxCordgOWaS425C6Qyg4JUBY4z+BPmHU422MaGr9xFjUiIodj1tAc7xxfjQkhUQq8gvEs8I/4A86i0u9j1kBNwqyjil0UeSvG5c+OwNeA2yW5D7G84CxALUtcfLu7AZ26tc87N3lLzILj0iabxTnApzGL0b6D6eX4irVQKEO5Z2Jc/m/ELFY7w9nvAxZiLIjewriHiLo4P4Mkru0xbmz6MR9iN2Yx19XEuA7HuHh5C+PL6ApiFog558cuBgwwfMHZGswCu3bMQr4atoQE6AKWAV9y0u7DltX43ThhEhJc4xuYRWD/xvja2iXmmGIWqb7OFg/QkqAc6zEL7cY52z/ErI6OLny9Hrgp9vqYnu9mjDPQ6AKzXZzyus8p/y7MQsGDkjy/pOEJXFwnq3AKcdffz8lPYvY9AVzn/P8T4NaYY7s4ZbpHkvtJGqqALSEhLnHu9x3gzJjjTc7z+wizIPU6ErjbT3DNxTjvfcy+e3EW8zrbLcC7zv+jcTyNxxz/NXBjkvz3wLi96WTLot1tY45vFcohST53Y7xAPIpZsHikU54PYr7T1cQshnXeo/ud59iFWXG/N8ZydgPmez8q7tk8jPkO3gC+4exPFlbCh3F/9I4j//VAdRLZr8bEmvmN83zOSXY+W0IFRN/Xu508DmVLXfUKEIjJP4gJ3/E359w9MY5tn3TuZwXOAvOYsrwVswC7C7Mge4+Y45+IOfc9nIXhePDNDCuXLCvt1ZiWz/cxLbd/AW9mk1ehfhjvtl2YlfGjMB/rdOfYQkyLtRHTslsJnO0cO8N58b7hvBznYVaGR40r7mZ4ZX41W1bT7+u8QDMxq31/4uSVibJpw6xGbnD2ncSWynMW5kPcOUbW+BXnQ9fABMH6ADgA43LlZuDZmLSKabVv65TX+8St6I5J+yxOfBhMhfsmW5Tws2xRgrHXH3Z/MeXVg3FmWY3xvvt8kmsejQkOtS1G8ewTc+/prnMxpkc+wbn3+cBvnWPfBP4P43uqGuMba2yC609la2XzJFtWtP8vcFvMsV2dMj0+yf18AVM5CyYUwccYJ5nRexjAuCga5ZTPx8B2zvHfYT7+0RgluD7+2Se5ZiJl8wpOLCZne5wjdxMwA9gcl/47mLVDifLfExOKoA7jYPRZtjQ8pmAq/ah3hkkkV8R3YxoUn8a869s4z/5KzLe0O0ZBHx33Hh2NaZQtxNRTlzvl9w1gdUz+z2A8ZNdjwj28j6P4SBBrBuNJe75T3jtgFPw3k8h+NeY7P8GRvSHV+Wz93e+KqcQ/75z/OWd7vHM8iPEq8gnnXn1OuZ7pbB+A+c4/EVOW/wYOdo4vAn7nHGvEKMBLnLJoBA7x6psZVi6ZVtrOhZ5LsK8um7wK9cO0cP6QYH81ptcQ60rim0DQ+f8M4I2YY9EAVjvFPMhkyubK6MOJObePzJTNWWnuqw2nMiO9srkL+HHMsTGYj2KSs63AzJjj9wHfT3Ld63Bc2GBccnwb45crvtcTe/1h9xdTXn+N2d6XuMot5thnMQ2BQ4lxjeLyOq8R04rGeMztd+Q/ixRB7mLOGYWp4L7r/H+U8zwfd463YD7yaZgKZj6mxXqqy3f0j8C3Y+5hM8Nd0Wxw7r3akf0/Yo7NiX/2Sa6RSNm8SUyjwrk3xSiDz+D0cmKOfwPn+3BxvROApc7/ezr3cCQwKs15dwMLY7YPAdYm+KZ/FfMePRlz7IuYhkG1s93o3NO2mMbbIMODFs5lS6/iaoaPTuyIqSMaYvadCjydRParGd6IS3l+/PuK6f39Oi7Px4HTnf+DwLUxx2YRVyc7795VMWX5y5hjnweWx8ixNMl95PzNxP6yHetrE5Fvx+5QEw+llJmI+ajiGYdpKcV6Y32LLa7OIUPX9jEMczfvnJupw8xhwZFE5DQRaZMtLur3w32ckV2IuU81noQ7SXKvmJZ0svt8BvORHIAZsngS0zo/FKOcP0hyXiLir1mfaJxbcwtPsBs5hlNQY+p/AqZH8i6mNXgfZrgLVW3FOJJ8EFPOazC96YTGCpI6VAGYcfxYzwPR5zGeLaGvo2TjTThKfCiL6P9dCY5Fj3eRABHZQUyIhvUi8hFmqGUcgGYYyoHh97cbsIsMD93xA4bHM3ov5v/NwAe6xZVW1BnrGMx3EI0nFCX+m49lN4wCfifm2vMxPRS3smdy/m7ASXH3OpPhIQXi8z8kLv3XGO6NPdl3naxejObrWQiSbJXNjhjPxm+LyCMicoOInJRlXoViHWbIIp4PMNp6t5h9zZhhiVwZ5m5eTDTNppjjblz1a8z5u2G8Bl8ANKlxb/9Ptri3163OHs7bxNynGN92TWR3r0swwyJfAp5R1WWYcvsCRhElIp18aVF34QkSXSfbcArx129X1cNVtUlVj8YM57wQc/xWVd1LjWeNBzFK4Z/x+Uj6UAWpeB8zxBYbvqDZxXnJGBbKwvn/PTVeplcCNTLcFdX+JI8IOxdT/tNUdSzwdWLuSZOHckhE7HNchxkGi31+jar6eXe3OIy3ge1luCf22G8+/v1Zh+mZjIu59lhVTWUkES97Juevw/RsYu91tA43worP/5m49GNU1Y3nimT1YvRYzt9MlKyUjaqerKr7YGJYXIl5IQ9OfVbRWQQcKSIni3HH3SQi052Wz30YV/uNToX+/zAtslx5ABMi4FMiEvXQGluZtAGfd0xLd8K0+lIxGvOSvQ8gImcyPM7Oe8AE51qJuBfjGn66U9nNAf6uqmsyuiuGemkvYzwgR5XLEswQZDJl8x7QJCK+TK8HGYUnSHSdbMMpxMswTYwr+W1E5DuY1ubdzrF6EdnPsRBqBu4Afq6qHybIKl2ogqQ47+xDwNWOHPuSwn2/I1utGPf3AoxyZI1+/wuBs0VkXzH+3H4YvSdV3eRc61ox4TI+jbG2+3WSSzXiGKiIyK7ENAYkdSiHdLwAfCQi3xORBjEhBfYTkU+6PH8INaGUlwBznXKYBpyNqSMgLqyEmkiyTwD/KyJjxfiI3MOxPnRzvUzP/w2m3jjauc96Meb8yWLXPALsLSL/JSbsxyjnW9nHhXiPADuJyMUiUufUgYc4xzz5ZqLkZDKnJrzwP1T1HlX1PACWl6jqWswwxSWYybI2trTmLsRUXqswY9r3Ags8uOa/nLx/h+nldGHGrKNDjr/GTM6uwbyMv0+T3zLMJHQI80FMxVikRHkK0+J8V0S2GsZyhnmuwLSo38G0aFKFJE7HM5jhgRdithsxk8KJ5F+OMShZ5XTNM12DMhbTs/uQLRZ+P3F5nZ9jrI+eEJEuzMRn9KPaCdMwiLr8f4bkjY3/wpTdBswczedihpDrMe9ON6ZMQpjy3gpnCOciTEPnQ0woiYfdFILDBZihkHcxiuFXadI/gancP4VRgpuBwxxZ/oIJI/A0plzfwgwHRvkWZg5qA6Zcz3Pe7URcgxlaDWOsnx6KOVaHmdf7wJF7B8xQWFocBftFzGT+aiePX2Imx7PhVMyc1NuYdYJXqeqTzrH7nb+dIhI1FT8N00BYhnleD+AuUmYU1+c7yvB4TNm8j+lhXEqS+tp5l47CfMtvY8r2R5jyTolz7ucwZfsuxhI1GjHUq28GyMFdjSVzxETu2wjspaqriyyOxWKxFAzPFgM5w0CWOETki85Qx2hMK/xVTE/GYrFYRgxerjy9y8O8KonjMV3bt4G9gFPUdictFssIIxevzzWxZpkispuq5mJ+abFYLJYKJeOejYh8Q0RWAOucydenRORQzCSbxWKxWCxbkZFzRxH5HnAQcLiqvuvsOwpjIZSLnb/FYrFYKphMI3UuB6ZqXNA0EdkR4xPqRI/l85xx48bppEmTii2GxWKxlBUvv/zyB6o6PtvzM3ZbH69onH3vicjN2QpRSCZNmsRLL71UbDEsFoulrBCRnObkM52zeVNEvpBAiGuB1lwEsVgsFkvlkmnP5lvAg2Ii/b2CWcH8ecxq/BWeSmaxWCyWiiGjno1j2vxJzJqaTZi1I19X1dOAe7wXz2IpDWbMmMHuu+9OKBQqtigWS1mSqTWaOAsS/+L8hlDVH8WlsVgqglAoRHt7O5FIhJaWFlpbW/H7/cUWy2IpKzKds3laRC50PNoO4XiU/ayI3EMa77MWS7kRDAaJRCIA9PX1EQwGiyuQxVKGZDpncwwmQttvRWQyxqlkPSZy4BPAz1S1zUsBS4ZQCIJBCATAtmpHFIFAgKqqKiKRCLW1tQQCgWKLZLGUHRkpG1XtwcTtvs2JKTIOE8J3Yx5kKx1CIWhpgb4+qK2F1larcEYQfr+fadOmEQ6HWbRokR1Cs1iyIGtHnE6ktncqXtGA6dH09cHgoPlrh1FGHD6fj+bmZqtoLGVJ6I5XmXt0kNAdrxZNhowXdY5IAgHTo4n2bLwYRrHDcmWFnaexpMPtJ+063R2vEnywk8CJTfhnT80pTcs396CPfah9oo9WXk2aNq+o6oj6HXjggZoVS5aozplj/ubKkiWqDQ2q1dXmb7I8vbymxWLJCjef4ZIlqg11A1otg9pQN5Dyk3aVbn67NrBJq+nXBjbpkvntWaVRVZ1z1NNaTb+CajV9Oueop9PccWKAlzSHujfrYTQnjvbI6Rn5/XDZZd70QtwMy0Xnia64wvy16zssFk8JhWDu3NSfVigELUcMcsXlEVqOGEyaNrjwLfp6lUGtoq83QnBhYs8urtM92EkftQxSQx+jCD7YmVUagMCJTdTSRzX91NJP4MSm5DecR7JSNiJyAfAe8JaIvCIi53grVoUTHZarrk4+LJfJPJGbr8aSkEAgYK3LKpB0n4TXSiTAM8MrdJ7JLZ0LBeFWifhnT6V1/ptcd9TfaJ3/ZnGG0Mh+zuYSjPfnd0VkZ2COiExQ1au9E62C8fuNRVuqgVu380SZWMrZeSJLmePmFY4qkr4+obZWaX26equ0RonsyiDV9PX2E1zYgd+/21Z5GeXwFfrQGOVw2lbp/KftReuCzxPs/zSBUX/Df9rchLK5Tjd7Kq2kno9xkyY2rX920sOFIZuxN+BVoCpmuwpoz2U8r1C/rOdsioGbweI5c8zcD5i/c+Ykz8vNPJHb61YIhx9+uB5++OEVf81Sx8t5kTnnrhk+R3Humq3zOvceZ76jz8x3nHtP0osuqT1c58gPdEnt4d58N2X6fZHjnE22PZvbgftF5Huq+gYmcNrHHug+Syx+f/oeiNseUKJhuUR52zVFlgLjpicC3vZG3PYw8PvxB+fiDwYhMDf1t+Dme80kXYWRlbJR1dtE5F3glyIyDRgLLBKRk4A2VX3dSyEtKXAzJAfeKyWoiGG5trabXKWLFpcXFtDhcJhwOEwoFKr4dTtuXpGiDGnlQ4lYUpNLt8j0rKgBpmF8ov0MeCrXPPP5K5VhtCVLluicOXN0SSG70q7HKlyaZbsdlvNKtjxwYPVzemD1c2nTHX64+eXKkiVLtKqqSgFtaGgo7PMvMK7NfIs5pGVxDTkOoxW98i/0rxSUzZIlS7ShoUGrq6tTVjhFUUiq2j5/vj591FHaPn9+8kSZzBV5peDywOflIf0iDyVdoxBl+nTVyZNzF23OnDkKKKDV1dU6J1m5aXHqS9fXdJHQzdxJNC+rREqfXJXNyFknU0IEg0E2b94MbPEiHD+cEgqFmDlzJpFIhIaGhoK5tQ+FQrRcfDF9fX3UPvccrVOnJr6um2E5t/M/mQzdeUjojld5Wo+mj1r++s3kK6tDIWhvh0jE3E4u01hunXoWY+rM9TVDIUKBy5xhqsvwBxMPQ7kd9ir3Ia3+/n46Ojro6ekptiieUF9fz4QJExg1apSn+VplUwQCgQANDQ2mQk9S4QSDQUQESK6QooRCIYLBIIFAIKVCcpPOjSIE3M0VBYPg5JVSiWTiDshDnyBmUdw+zqI4JfhgZ0Lz0GDQKJp0t+EGt049i6F/g0HMmpKI0NerBIOS2IZk4eu09D1KH7XU9vXRuvCBhPfhehIeSlKJuKWjo4PGxkYmTZo09M2WK6pKZ2cnHR0dTJ482fvMR9KvFIbRVNMPkWUy1OY2nZu5Arf5ucJrc+tM5pOqqswQX4p0W9x99KV09+EyO9e4MX3Ox8hiuiJ2Wx6uh8fcXLQCWLZsmUYikWKL4RmRSESXLVu21X5G8jCaiOwOXA74VPUrxZYnE/x+f8peiN/vp7W11bOeSKIAYAlboy6v6wq3lnLRtOmuFQwS6j2AYOQzBHqfM8MuyXpULroi/tlTOeLcP1CtcNn8PZMuivP7Ydo0CIdh0aLCNMAzKTo3uBki83c+QmvVn035Vj2Hv/MLwNZlEjhtN2p/NUhf3yC1tVUETtvacmzYjZRpjyUTMu3RREczStHBa756ZyWnbERkAXAcsEFV94vZfwzwc0ygtl+q6o2qugo4W0QeyLdcM2bMKHg8k3QKCdwNyWWSzu113RK47DIAgs7fXAg1HUdL5Ntm+CbSR2vTmySUMhAgVHu4Gb6p/hv+FPf6XtV4gLQuPHw+80s9cuduONM9ISAIBCDxnbrG1bBcIIC/7jr8fc87w5n/kzAvvx9an64ud6t3S4HJStk4vtEWqeqHHssDcDdwC7Aw5nrVwK3A54AO4EUReVhVl+Xh+luRjxj0Xi1R8fv93HTT33nwwU5OPLEJvz9Z69zDHksecLUeo3MqPSiK0FdVTbBzasIqOISfw/pbGVShXpRWqnOsqtPjtUGH1/kFAlBbM0hfBGprIBCo3jpRBt2pEdJhsXhItj2bnTAV/j+ABcDjzphezqjqsyIyKW73wcAbTk8GEfkdcDzgStmIyGxgNkBzc3PGMrkdgnKLl5ZGoRBcfPFU+vrguedg6tRUefnJtYUci1dDAW7LIxCA+gZx0knKdakDairTvoFCTa57+454nZ+fEK193yGohxHQEH7mkvBdsFqkIIykhb1RsvL6rKo/BPYC7gLOAF4XkTkisoeHssWyK7AuZrsD2FVEmkTkF8AMEUk6TqOqd6jqQap60Pjx4zO+eHQIqrq62pMY9G4dOrtx5pxJXm4jFrh1Ih0O78vatV9Nme7tt3ejvf0L3JEiQqDbe4g2vK+7LrWCNo8nAmha47bl9W/yz5rfE8rRY3Ym70gwGEyroKMm0oAn7xzBIP6qv3MZN+IfXGyjzRaR6EjJ6tWraWlpyfndi+W2225jv/32Y7fdduPmm2/2LF9PyMW6ANgfuAlYjvGXthT4cS55OvlOAv4Zs30SZp4muv1fwM3Z5J2tNZqXCyzdWBq5tYBya7Xktb/OJUtUoV9hQOuSrA6fP79dYZOTbpPOT2Ht5bXl1fTp5+nkyfPTGLctUahTqPJkNb/Xi3A9za+IC2crnUSWW6nIZGFvJjzwwAN6wgknaF9fn7799tva1NSk/f39WeVVMtZoInIRxj3NB8AvgUtVtV9EqoDXge9mpfmS0wFMjNmeALzt8TVS4uWkudslKlGjkFTrLIrlGm3hwrcwHc5qenv7WZjAn9WDD3YC+2BGa5UHH+xkdoJ1LF5bXgH4fMvw+ZbhT+FX3fQu+oGIJ0NVXg9T+r/1LfzhsDdhyPNRyJascLuwN1PmzZvHnXfeyahRo9h5550ZNWrU0FBsSZCNhgKuBXZLcmyfXLSfJu7Z1ACrgMlALfAK8Ils8s6+Z1PY5QLFWGeRyXXPPfcep9fSp7BJz03gz2pLz6YvZc9G1b3bfS/TeemnzPPn5fXiHkveyLRno6o6ffp0nTx5sme94L6+Pt12222Htt9++23db7/9dMOGDXrGGWfounXr9Mwzz9xqu6+vL2F++ejZZBsWuk5Vh4WsE5EfOcrrtSzzjObzW4zN5xQR6RCRs1V1ALgAeBx4DbhPVf+Vy3UyoRgRmt3OT2SaZ7rI1m6ve9ppeyHyOeBKams/z2mn7bVVmtmzp7LXXuex3XY3MX/+m8wuUoTAZERX80+ePDlna69MAqu6zjB+rZClYvD5fDQ3N3s2WrJs2TLC4TCrVq0iEolw2WWXceGFFzJ+/Hiam5u55JJLmDdv3lbbXrukSUW21mifA74Xt+/YBPsyRlVPTbL/UeDRXPPPBrdeV7ymWIZBbq7r9/vZf/+PCYd/n3Lt0S67vMUuu7zF7NmXeiKb11Y8Pp8Pn8+Xc16ZeNxxnWFDg4cZWiqZpUuX8rWvfY1TTz2VTZs28eUvf5nZs2fT3d3NqlWrqKmpYcyYMVttF5KMlI2InAd8C9hdRNpjDjUCf/NSsFLCfveJ8aqiBndKJB/rndqmt+V0fpRMpkRcxcYp8TkW7xewjiy89hzQ1tbGcccdx6xZs4b2DQwMcNFFF3H99ddz33338de//pV77713aDv6/ApFpj2be4HHgLnA92P2d6nqvz2TqsQo8e++7HGrRLxee+I1nvdES3TNS7E8kluS09bWxnnnnTdsX01NDQsWLADg0kvNyMKRRx45bLuQZKRsVDUMhIGEQ12VTIl+9xWBWyWSLyseV7hxcVABkUvdUOpKfyRSij7W4sl0GG2xqs4UkS6MnfjQIUBVdayn0llGBG6ViFv3/J4TCsHMmWbCvqEhseWEmzQO4bD5hUKF1UleeXzIxM+exRIl057NTOdvY37EsZQTXrWmMlEiXs4TucaNF2mXnqa9DMRWLErdz56lNMnK9FlE/ltEdvVaGMvIxWtTULcMVA/QU99DaF0Ke/ZAABzXMUktRKJWJNXVKa1IimnRHA6HWbt2rSfuUfx+P5dddplVNBbXZLvOZizwhIg8JyLni8iOXgplsRSC0LoQmxo30VvfS8vCluQKx++HxYthzpzkXRGXC5Tc6K18kE9/XBaLG7J1xHmNqn4COB/YBXhGRP7qqWQWSw64acUH1wTNPwJ9g31bthPhdkVsmjTRQGyTJxd2CC3RpH4iAoGAnYMpAIG7AwTuDhRbjIKSbc8mygbgXaAT2CF3cSyW3HHbig9MCph/FGqra7ds5xmfD5qbC28c4KkXaYslQ7KdszlPRILAU8A44BuqOs1LwSyVgRt3+vm4pptWvH+in9Fdo6nrqaP1tFb8Eyt3/sGtax4v53Uslliy7dk0A99W1X1V9SotUMRMi8UNmbTiawZrqO+pr2hFEyWdEYad1ykc4d4wa8NrUxumZEEpx7PJSNmIyGLn3wuB50TkI+fXJSIfeS+exZI5XjrYHEm47RFaciO0LkT7e+2s3rg6tWFKhjz44IM8+eSTLF26lOeff55rrrmGgYEBT/L2gmzX2RTWg5vF4uC2AizKepwyp6geGkYQwTVBIuoodccwxYuedanHs8l2zuZHbvZZLJbywfYIC0NgUoAqcYZ5PTJM6e/vp729nb333huAd955h3HjxrFo0SIee+wxVJWzzjqL22+/fdj25qg7+wJQciEGLBZLfnDTK7Q9wvzjn+hn2o7TCPeEWfTlRZ70amLj2UyaNGkons1hhx3GggULWL9+PbNmzWLPPfcctt3Q0ODBHbnDhhiwWAqM26mQ0LoQwTVBApMCZWvAMGPGjML7sisDfHU+fHU+z55rsng20WPhcJhzzjkn4XahsCEGSphKqGws2RFaF2Lmr2YS0QgNNQ1laZqdj/hDlsQkimcTpaamhiuvvDLpdqGwIQZKlNC6EC0LW+gb7KO2urYsK5ti4lVQtGKRr0nktNf10ALNhiJITvCMoKf5JYpnEw6Hufzyyzn99NPZYYcdttouNFnN2YjIScBfVLVLRH4IHABcr6r/8FS6CmbG/Bkpx2yDa4JsHjCTd4WsbCylQWBSgIaahqHGRqG8G3iJtW4rHIkaCT6fj1tuuSXpdqHJdlHnFY6imQkcDdwD3O6dWKVHaF2Iuc/N9cQm3o2dfbSyqZbqsq1sLNnjn+in9bRWrjviurLt1VrrNkss2VqjDTp/vwDcrqp/EpGrvRGp9PB6/NzNEEm0srFzNqVB1GmiJ8Mf0RZ+miEr/0R/2T/3pUuXFlsES4mQbc9mvYjcAcwCHhWRuhzyKnkSKYdccGtn75/o57LPXFb2FU4lkC/3IhbDjBkz2H333a2LnBi6u7t555136O7uLrYonpCtgjgJY5V2lKpuBLYDvuOVUKWG14uw/BP9LD5zMXM+O6egQyReDgWOJLx2LxIaG2Zus1VcUTLxyRYKhZg7d27FK6Xu7m5WrlzJ+vXrWblyZUUonEzX2Sx2XNZsANTZByDO9livBSwFosrByyGtQg+RFMuU1u3wUyn74fLSMiy0LsQR+7fRXwV1C1vKdj7GS9xarYVCIWbOnEkkEqGhoSHlPFAoFCrrsNVdXV1DZRKJROjq6mLMmPL2Epatb7TG/IhTupT7+HmxTGkrgWjPNqKRnHu2wTVB+qsgUmWfQxS3VmteKyW3Csltuu7ubrq6umhsbEypGNyka2xsHCqTqqoqGhvLv8rN1kDAUmZ4WWFGGSmLTr10LxKYFKBOq+iLVqzWynDIai2dpwEvlVIoFKKlpYW+vj5qa2tTKiQ36aLDXlHlsPfeeydUJNF0e33jG6gI3cFgwnRjxoxh7733dqW8ygUv1tlcAcygCOtsRGR34HLAp6pfKeS1yw2vhwLdDsuFe8OEe8KE1oXKWiF55V7EP9FPa9s0gtuGCdzgjV+sSsCN1VomSqmhoWFIQSRSSsFgcMgJZapeUjAYpK+vj8HBwZTp3A57xaZDNeXw2JgxYypCyUTJtmdzhare76yzOQr4CWadzSFuMxCRBcBxwAZV3S9m/zHAz4Fq4JeqemOyPFR1FXC2iDyQ3W2MLLwcCnQzLBedWI9ohBY7PzGEfx34/wl0ABOLLU154VYptba2phz6cqOQoulqa2vTpnM77BVNB4BIRQyPuaWY62zuBm4BFkZ3iEg1cCvGq3QH8KKIPIxRPHPjzj9LVTdkLrrFC9wMy1XSPJFn7kVCIWhvh0gEWlqgtRXKcALb03VHecDv96ecX3GjkDJJ53bYK5qutq+Pmq4uql591dPnf9ttt3HbbbfR1dXFd77zHS688ELP8s6VbJXNehGZj1EKP8pmnY2qPisik+J2Hwy84fRYEJHfAcer6lxML8hSIriZx8hknqjUKy/PCAaNogHo6zPbZahsKmF4NJ1CyjSd22GvMa++CsuWed7giI3U+cEHHzB16lTOO+88ampKY2o+23U2JwOPA0c762y2By71QJ5dgXUx2x3OvoSISJOI/AKYISKXpUg3W0ReEpGX3n//fQ/EtAAs/eZSVn17VdLKJqqQJm872Q6hRQkEIDqMUlu7xZtAIS59d2BIqedCvsIajxgSNTg8YN68efzoRz+qrEidwGZgNFu8P48CNnogjyTYp8kSq2qnqp6rqns4vZ9k6e5Q1YNU9aDx48d7IGZpMWP+DHb/+e4l+dH76nw0+5qtooni98O0aTB5csGH0LzyguC1R40RRx4aHMkidYbDYc4880xWrFjBKaecws0331y0SJ3ZKpvbgEPZomy6MHMtuRI/ZToBeNuDfCuWSmlljih3MD4fNDcXVNF4+Z7kI6zxiMLDBkfUpc1LL700FKkzEokMReocP348O+20Exd++0K+dem3aN6jmaeeeoq77rqr4JE6s1U2h6jq+UAPgKp+CNR6IM+LwF4iMllEaoFTgIc9yLdiqYRWZqUozFLGy/fEDo96gIsGR3dfN+90vUN3X2JXNbEubR5//HFmzZrFqaeeyrRp02hubmb27Nl0d3ez4vUVVDVUsc1O27Dr/rvy0j9eoq2tjaOPPjpfd5eQbJVNv2M5FnVZMx7IaHBQRH4LhIApItIhImer6gBwAWY+6DXgPlX9V5YyjggqoZVZCQqz1PH6PbHDo/mlu6+blZ0rWd+1npWdKxMqnNg1OytWrODwww/n73//O//85z+59tprGRgY4KKLLuKSH17ClE9M4eXQywDUjCqDSJ0xzAP+AOwgIjcAXwF+mEkGqpow2qeqPgo8mqVcI458+G0rNNZqLf946QUBbPnnwooPVsAD85kybkrSNF29XUMNsIhG6OrtYkztcEu3xsZGpE7QUcrKN1Zy6aXDbbRqampYsGAB3X3d1H+rnq6PuvjxD3/MGWeeUR6ROsV43nwWeBlowUzqn6Cqr3ksm8Ul5e63zeuKsFiUuiL0yguCJTcGdZDByCDdfd1bKZAojXWNVHWbBliVVNFYl2DxZy3QZP694093sPe4vRPmNaZ2DFPGTaGrsYsF8xckvWa+yVjZqKqKyB9V9UBgeR5kshSZdCGr84HXFeFA9QADowbKeh2IpfLo7uvm4/6PAVjZuZK9m/ZOWPmPqR3D3k1709XbRWNdY8I0Xb1dqGOsq2jC3k9sfsVSMlGynbN5XkQ+6akklpKgEibrQ+tCbGrcRG99b9neQz4InhEs2V7XSKGrt2vo/+jwWDLG1I5h58adU/d+nHm4pL2fEiJbZXMEEBKRN0WkXUReFZF2LwWzFIdKmKwfklnK9x4spYOXQQdjFUKuCiLa+9m1cdekPaRSIlsDgWM9lcJSMuQjFEGhGZJZobamPO9hpFGq811eBx0cUzuGbUZtw2BkkMnbTc5ZQZTC8JhbslI2qvqW14JYSoNKmKz3T/Qzums0A6MGaD3brgMpJl4rEbf5eXXdfDiTrZZqqqury0ZJeEVpeGizlBRLv5nejXupUzNYQ81gjVU0lpwITArQUNNA32Bf2fb0SwWrbCwWyxClOpxVLPwT/bSe1urpOrYp46YMuUPzyAdnWZBtpM4fqer30u2zWCyWcqfc17GVCtlao30uwT5rNGAZ0YwoZ6IWS4ZkpGxE5DwReRXjz6w95rcaeDU/IlospU8lrE8aaXgV3ycbwmFYu9YEbvWS2267jf3224/ddtuNm2++2dvMcyTTYbR7gccwIZq/H7O/S1X/7ZlUFkuZUUkhsC35JV+RwUs9UmdGUqhqGAiLyJnAl4FJ0TxEBFW91nMJLSOCcp+QroT1SZbCkK/I4PPmzePOO++suEidfwSOBwaATTE/i2VEYmO8JMbOY21NPiKDJ4vUuWjRomGROdetW8eZZ55JR0cHZ511Fv39/blf3CXZ9q8mqOoxnkpisZQ51qvycKLzWBGN0LKwpeKV8IoPVgCkDB0AWwJ1hsOwaJE3vZply5YNReqcNGnSUKTOww47jAULFrB+/XpmzZrFxIkTaW5u5pJLLuGuu+5i1KhRuV/cJdn2bJaIyFRPJbFYLBVFJfjZyxdeRwZfunQpX/va17aK1LnHHnuwdOnSocic3d3drFq1ipqaGsaMKawHg2x7NjOBM0VkFdCLiWmjqjrNM8ksFktZU8x5rHBvmHBPeMSEmGhra+O4445j1qxZWx2rqTGROaORO6+//nruu+8+gsEgAS/G8FxiHXFaLJa8UCw/e+UwfOe154C2tjbOO++8YfvC4TCXX345p59++lBkzgULFgBsFdWzEGSrbNYCXwN2V9VrRaQZ2AmwDjotFssQxZjHcmuGXkm9n2AC7eXz+bjlllsKL0wSsp2zuQ3wA6c6213ArZ5IZLFYRhxeWq1Fh++ApMN3dhFu4clW2RyiqucDPQCq+iEmIrbFYrFkRCYVvxul5MYMPR/GC4M6SN9gH9193TnnVYlkq2z6RaQaTABsERkPlM7qIUvFYtdtVB5uK/5MlJKvzkezrznp8Jib3k8mdPd183H/x/QO9rKyc6VVOAnIVtnMA/4A7CAiNwCLgTmeSWWxJMAOfVQmbit+L3sjXi/C/ajno6H/Ixqhq7crp/yKiarmJd+slI2qLgK+i/GR9g5wgqre76VgFks8dt1G/ilGz9Ftxe91byRd78ct9fX1RD6OOOM8UCVVNNY15pRnsVBVOjs7qa+v9zzvjK3RREQwHgSWA8s9l8hiSYL1P5Zfimky7MZqrVRDlk+YMIGOjg42vr8RFMZtM451G9cVW6ysqa+vZ8KECZ7nm7GyUVUVkT8CB3oujcWSglKtbCqFcvBcXYougUaNGsXkyZM585kzgfJ3Kpsvsp2zeV5EPumpJBaLC7wa+rBsjdfDVBZLLNku6jwC+KaIvIXx9mzd1VgsZY7tOVryiXVXY7FYhijFYSpLZSDZmrmJyHbAXsCQ2YKqPuuRXHlDRMLA6zG7fEDY5fY44AMPxYm/lhfpk6VJtN/NvtjtfJZFMnlySZ/qeDblMZLejUT7K+lb8frdiN8up7JIlSZ2/26qOj6D6w5HVTP+AecArwIfAk8Dm4Gnssmr0D/gjmy3gZfyKYsX6ZOlSbTfzb64+89bWeSjPFIdz6Y8RtK7ke5+y708vH43Ur0rpV4W2bwb2fyyNRD4NvBJ4C1VPQKYAbyfZV6F5v9y3M6nLF6kT5Ym0X43+/4vxTGv8bo8Uh3PpjxG0ruRaH8lfStevxvx2+VUFqnSeHYfWQ2jiciLqvpJEWnD+EnrFZE2VZ3ulWCliIi8pKoHFVuOUsCWxXBseQzHlscWbFkYsjUQ6BCRbYE/Ak+KyIfA214JVcLcUWwBSghbFsOx5TEcWx5bsGVBDgYCQxmIHI6ZRPqLqvZ5IpXFYrFYKoqsejYiUg98CxMeWjGOOLOd/7FYLBZLhZPtnM19mIBpv3F2nQpsp6oneSibxWKxWCqEbJXNK6q6f7p9FovFYrFA9kNfS0Xk0OiGiBwC/M0bkSwWi8VSaWTbs3kNmAKsdXY1A69honWqjkAfaSKyD2b90TigVVVvL7JIRUNETgC+AOwA3KqqTxRXouIiIrsDlwM+Vf1KseUpNCIyGrgN6AOCauJhjVhG6vuQbc/mGGAycLjzmwx8HjgO+KI3ohUOEVkgIhtE5J9x+48RkRUi8oaIfD9VHqr6mqqeC5wMlK1NvUdl8UdV/QZwBjArj+LmHY/KY5Wqnp1fSQtLhuXyZeAB5534z4ILWwAyKY9KfB/ckK2yeQ84EfgZ8FPMy/Seqr6lqm95JVwBuRujQIcQkWrgVozT0X2BU0VkXxGZKiKPxP12cM75T4xlXmthxfeUu/GgLBx+6JxXztyNd+VRSdyNy3IBJgDRaGKDBZSxkNyN+/IYkWS7qHMhxhrtZmf7VODXQFlao6nqsyIyKW73wcAbqroKQER+BxyvqnMxPbhE+TwMPCwifwbuzaPIecOLsnCiud4IPKaq/8izyHnFq3ej0sikXIAOjMJpo0KXSGRYHssKLF5JkO2Dn6KqZ6vq085vNrC3l4KVALuypTUG5oPZNVliEQmIyDwRmQ88mm/hCkxGZQFcCBwJfEVEzs2nYEUi03ejSUR+AcwQkcvyLVwRSVYuDwEnisjt5N+/XimRsDxG0PswjGx7NktF5FBVfR4q1hpNEuxLak2hqkEgmC9hikymZTEPmJc/cYpOpuXRCVSi0o0nYbmo6ibgzEILUwIkK4+R8j4MI1tlcwhwmogMs0YTkVepHGu0DmBizPYERob/t0TYshiOLY/E2HIZji2PGLJVNsekT1L2vAjsJSKTgfXAKcBXiytS0bBlMRxbHomx5TIcWx4xZDVn41icfQTsCOwW/ZWrNZqI/BYIAVNEpENEzlbVAeAC4HHMGqL7VPVfxZSzENiyGI4tj8TYchmOLY/0ZLuo8xzMAsaohcmhQEhVP+updBaLxWKpCLIdRotG6nxeVY8Qkf8ArvFOrPwxbtw4nTRpUlFl2NS3ieWdy0FBRJjSNIXRtaOzzm9F5woApjRN8UQ+t/ktW7aMwcFBdt99d0aPTix/27o2AKZPnJ76miuca07x5h7cXPfd7ndZ/9F6syGwa+Ou7DRmJ0+unwo39+r1O1Is3L5LxUhXrO+mXHn55Zc/UNXx2Z6frbLpUdUeEUFE6lR1uYiURQlPmjSJl156qagyzH1uLj946gcAVEkVpx1xGpd9JnsLyBnzZxDuCXPzl2/GP9GfNF3g7gAAwTOCKfNzky4UCjFz5kwikQhvvfUWra2t+P1bX3vbi7cF4KWbUpd5IOBcM5haNre4uW5oXYhP3fUpABpGNXD/afenLD+vcHOvXr8jxcLLd87rdG7zcovX+ZUaIpLTFEm262ziI3X+iRFsZZEpgUkBqsQUfW11LYFJgazzCq0L0f5eO6s3rqZlYQuhdSGPpExNMBgkEokA0NfX55mS8IqB6gF66ntSlod/op/RXaOp66mj9bTWgigat3j5jlgspUC2BgJfUtWNqno1cAVwF3CCh3JVNP6JfqbtOI3J207OuZILrgkSUafSH+wjuCbokZSpCQQCVFU5lWFt7VBrPVvC4TBr164lFMpdWYbWhdjUuIne+t60CrhmsIb6nvqSUjTg7TtisZQCObuOUNVnVPVhGxI6M3x1Ppp9zTlXIsVqAfv9fqZNm8bkyZOTDqEBDLwzQM+LPSmVSCgUor29ndWrV9PS0pKzwhlSuFJYBew1Xr0jFkspkO2cTUXR399PR0cHPT09BbvmVZ+4CoDXXnstp3y2ZVuePOZJIhph3DbjqOuuY/Xq1UyYMIFRo0Z5IWpSfD4fPp8vqaIJhUJsemgTDEJLS0tSpRQMxvTOnCG5ZHm6ITApYNbzK9TU1NghKIulBMhY2TixKXpUtWK8t3Z0dNDY2MikSZMwPiTzT9UHpjcyZVzudhWxeakqnZ2ddHR0MHny5GHpwr1hwj1hQutCBWktB4NB4+NXUyuRQCAA1cCgN0NydAD3AM2gb6vx0jYxzTkFpG16W7FFsFgKTtphNBGpEpGvisifRWQDsBx4R0T+JSL/IyJ75V/M/NLT00NTU1PBFE0+ERGampq26qXlw5CgbXpbyopzSIlIaiXi9/sZ/eXR1PnrUg7JuSUYDJqwfs/B4JrBkjNecEvwjGDFWjZZRh5u5myeBvYALgN2UtWJqroD8BngeeBGEfl6HmUsCJWgaKIkupdMDAnCvWHWhtfmrJAyUSI1O9dQ/8n6nBUNuFdyFoulcLhRNkeq6nWq2q7q1FaAqv5bVR9U1ROB3+dPxMpkyrgpngyhucWtIYHXPSAvlYhbvO4peU3vxl66P+jmjsfuKLYoZYdXDSFL4UmrbFS1H0BEThKRRuf/K0TkIRE5IDaNpXRxa0pbLFNqr/FayQXuDgwt2suFOx67g54JPQzuMcg3//ZNq3AyoFhryizekInp8xWq2iUiM4GjMFOwt+dHrNInEAiUzPCM216SG1Nau5gwvzz48oNmiK/K/B58+cFii1Q2VEpDaKSSibKJWp99AbhdVf8E1HovkuWhhx7ioosuKsq17WLC/HLigSeaL2kQiDjbFlfDY5k0hEp5uM2rXnK5kYmyWe+EPJ4FPCoidRmeb3HJ0qVLOeCAA4p2fbeLCd24hPGaYn2obTfeRNuNN+Wcz+xjZ1PfUU/1qmrmf3o+s4+dnbtwZY7b4TG3DSG3+blVSCNVOXhNJsriZExchqNVdSOwPXBpPoQqB7x0rxJl5cqVfO5zn+Omm27i6quv5qabbvIsb/DWlDYTlzCW4dRtW8eYcWOsonHIZHjMTUPITX52/qfwZKJsNgOjgVOd7VHARq8FKge8dq8C0Nvby8knn8xPf/pTxo8fz/PPP8+1115bUK8GmVApLmEsxcfreUI3+dn5n8KTibK5DRMkLapsuoBbPZeoDMiHx+Mnn3yS/fffn1122YWxY8ey0047UV9fz+BgaTpqGPqA1RoSZMr0tulMb5tebDFKBq/nCd3kVynzP+VEJsrmEFU9H+gBUNUPGaEGAl57PAZoa2tj6tSpvPLKK0ybNo0NGzbQ2NjIxx9/zJlnnklHRwdnnXUW77///rDt/v7iWJ37J/qZvvN0Jm+XuoIoxryOpfzw2ulouvy8nv+xpCcTZdMvItUYF4eIyHggkvqUysStx+NMGDt2LMuXL6e9vZ1p06Zx1VVXcf755zN+/Hiam5u55JJLmDdv3lbb+Xa2mYp0H7Sd17GUMl7N/1jckYmymQf8AdhBRG4AFgNzvBRGRBaIyAYR+WfMvu1F5EkRed35u13MsctE5A0RWSEiR3spSzp8Ph/Nzc2eLRr8+te/zuuvv861117L7bffzvbbb8+FF15Id3c3q1atoqamhjFjxmy1nQ+8MiSw8zqWcseuO/MO18pGVRcB3wXmAu8AJ6jq/R7LczdwTNy+7wOtqroX0OpsIyL7AqcAn3DOuc3peZUl22+/Pc888wwTJ07khRde4IYbbmBwcJCLLrqI66+/nunTp/PXv/512HapO5ispHmdgaZX6Nn3Lts7G2EUc91ZaF2Iuc/NrZh3LqMQA6q6HOP1OS+o6rMiMilu9/FAwPn/HiAIfM/Z/ztV7QVWi8gbwMFAQZ5MPir63t5eurq6aGpqAkwslgULFgBw6aXGyvzII48ctl3KRMMuD4waoPXs1B/qQPUAA6MGChb+IBNC60Js+uw3oaqfloV3VfRi1+j8Y6k3ZAqJr86Hr87n2TN3E+ojtC7EzF/NJKIRGmoakhs6OOt/ysE7uOuejRi+LiJXOtvNInJw/kQbYkdVfQfA+buDs39XYF1Mug5n31aIyGwReUlEXnr//ffzKmwu1NXVsWrVqmKL4Sluwi6X+txOcE0QqvqhKuLJcGBb2020td3khWgZUUoulkYqbg0O3M4VlZOlXKamz35Kx/Q5UUwATZRQVe9Q1YNU9aDx48fnWayRw0iZ2wlMCkBkFESqy3440FJc3CqRwKQADTUNVEvyd67cLOUyGUY7RFUPEJGlYEyfRaQQps/vicjOqvqOiOwMbHD2dzA8/uIE4O0CyGPJgOnTp6dNM2xup6b0KnP/RD+jn5rPwI4v0frzr1bsEJolMV4OUUUNDiIaSdlw8U/003paK8E1QQKTAgnfuUSKq5TfzUyUTbFMnx8GTgdudP7+KWb/vSLyU2AXYC/ghQLIY/GY6JqdcE+YRV9eVJofTIcfVn4OOnYpqRDTlvIianDg5l33T/S78tCeTnGVCpkom3jT568AP/RSGBH5LcYYYJyIdABXYZTMfSJyNibY70kAqvovEbkPWAYMAOeramkut7ekxetJWC8JhWDTuj1BhZYWaG2FXCze3fT2RhptbW3FFiFn3Ez8g3fveiaKqxRwpWzExBl+FngZaMHMl5ygqq95KYyqnprkUEuS9DcAN3gpg8USTzAIqABCX5/ZLrHgn65om95WbBFyplSVUnT+JKIRWha2FMxicek3l7pKVwpWa66UjaqqiPxRVQ8kj6bPpcyKFebvlMJFcrYUgO4HHjH/3JQ8TSAAiJo5pVqhkg26vFZIpaocvKbc5k+KQSbWaM+LyCfzJkmZYWNcjBz8fhg98Q3qxr2T8xBaqTPwzifpefG/SefIPBSCuXNJm85tfuWO9TSQnkyUzREYhfOmiLSLyKsi0p4vwUqNwUHo64Pubm/yW7FiS28pnmJG6rQkpqbhY+qbNpS1oklX8YdCsOmhP9EbupyWluSKJBSClha44grSpnOTnxvZMk1XaEo9wm0prMfJRNkcC+wOfBb4InCc87fi6e6Gjz+G3l5YudI7hZOMYkfqtJQGXi7CDIVg0/1/oXfJlUkr/mAQGGgArRmam0pEMAibN+tQAyxVOjf5ZaLk3Cqv8Bv7svaRr+acJhO89lztFaWyHicTa7REwdLDIvKyqrZ5JE9J0tW15f9IxGy7tTxJxuCg+XV3Q9Sf5sqVKzn//PN5/vnnaWpqYuPGjVx88cVp8+ruNjI1Nm7JqxBE68FUnk2GQimfkV9Z4nEzF+M1pezqxRg5VJHKyCEQAKo3w+AoamqqCAQSuxoMBKCqtpfIQA21tTVJ57ACAaDG5FdbOyppumAQGKwdppQS9SCHlFcaQ41QCNrm3AxaRcujia0HQyFo/5+fEukflTRNbNpg0NxPqp5t+I19CS+fTmhKaQ21lsp8UiY9m4OAczEuYXYFZmPMlO8Uke96L1rp0Ni45f+qKljWlVtLIVFPKdtInd3dsHw5rF/vXa8rEKCiJ8Fj0d5/Eun6ddpoq9N3ms70naYXRqg8MFTxSz+1tcmebwhj+Hklqi0kczPo98O0S/8fk7+0IGUl7fdD/VeOo/oz13LTva8mTRcIAFV9QD81NYOplVf15rTpjK5XYpVSojSR/lGg1Wl7XTNnwg9+kH7IsP1/fsrqh872pNcVzTPdvJibNIFJAaTjU/Dc96leP7No80mZKJsm4ABVvURVL8Eon/HAYRS83VpYxoyBbbaBujrYe2944b3cYlwk6iklitRZU1PPK68MplQgifJKhlslEg7D2rXpJ3/dpnND2403bekFeZDODaFQiEjnxehHv/QsvHep4vdD/WHHUT3xWm66KXHFHwwGIfI8cCODg4tT9tB8ey6j+bh7U7bgQ+tC9Oz9NIOB67no1U+maJS5U3Ju0xml1JdSsQYCUDWqH6oGUihfo4TEcYyVbsjQrfJqm3MLqx/8Rs7zYm7nzujwwz2t8NR1yMJWs10EMlE2zUBfzHY/sJuqbgZ6PZWqBKmuhtpao3hytTyJ7yk1Nm4dqXP16g00NDTy5z/fx69+9RhdXcpZZ53F7bffzmOPPYaq2a6p2bxVXslwoxxCIWhvh9WrXbTkXKQb2LwNPZ07pJ/49TidG4LBIGg/EPEsvLeXhMNh1q5d64kSDIVC9Dz7NIPrrufiiw9JmGcgEIBqQLyJQLvw2YXmnyro7e/dsh2HWyXnNp3fD6O/fDx1/jlJe15+PzSc8EXqDr0hZe8sEDDfffT7T9Xrcqu8qquGD2cmS9fXR8p5MTdpoul0YBRoDQP91SmHvfNJJsrmXow12lUichXwN+C3IjIas4p/xJCr5Ul8T2nMmK0jdV599VWcdNL5HHDAYSxdupj58+9i1qxZHHXUUSxevJi77jLb48c38B//AbvuuiWvRLhVDsGg6SGBi5ZcmnTRlfe9H+ycfuLXw3RuCQQCIKOAKk8qV7fKwZ3SD9He3s7q1as96XUFg0HYGfg09I7vTVhR+/1+Rn95NHX+Om8i0K4BBp1fxNlOgFsll5kyDGFCbyUvt5qdX6T+kz9L2Tvz+818znXXpZ7XcTu0mInySpfOa0WYd1TV9Q84EPg2cDFwUCbnlsrvwAMP1HiWLVu21b50HP6rw/XwXx2e8XlRli83vyidnZ162GGH6Xbbbae77767fuc7P9AXXojoiy+qfvrTx+rs2ecPpT322GP1/PPPT5DrFuLvac4cVTC/6mqznYglS1Srqky6hgaznW06c81I2mt6nU5VtWpcm8rY1Unl35LudpWxs3VJmoSHH25+yViyZIlWVVUpoA0NDUnzc1u+c+bMUczEg1ZXV+ucVDfrgvmPzlcuR7kS5XJ0/qPzE6bzfdunvm/70ubn5v1fsmSJshvKZ9DaPWpTlvHok0dr3afq0j4HN+mWLFmi1KBI6mfh9ppucVsnLFli3t10l3WTzsu80gG8pDnUva6t0RyXNfsAPlW9NhrPRlWt88ssiPdEEI3Uuf/++/PUU0/R1NQ0ZGU2dmwN11135VDampoarrzySjIhEDDDbJFI6laQ3w+LF6e3vvH7Ydo000pftChxOrcr771OFwpBpHM/Y42UxpeZ1O2H1O2Xcys+GAwScbp60SG5RHkm6hEmLrsAVVVVRCIRT3pdnWM6je2pQFVVldnOATeeAfx+P6MPGc1AxwBP//rplGVcs3MNNTvXpH0ObtIFg0HTm9LkzyIUCrHpoU0wCC0tLSl7cqFQiGAwSCAQ8CQMvN/vzlrNTTov88o3mZg+34bpDH8WuBYTz+ZBYER6FciHj6H4SJ2Dg2FuuOFyzj77dHbYYQfC4TCXX345p59utjPBrRKJpnXzYvp85pdKIY2e+AYDH4+h9eFdCpbOjZmv17hVDu6Vvp9p06YRDodZtGhRzpVcbBiHulF1OVskuY2s6laJeMnQcNtg8uE2NwoJjKKZOXMmkUiEhoYGb4YXRyjlEM9mxBAfqdPn83HLLbck3c4Ur1s3biYaaxo+pqbhY/z+XQqWzvSAIqBCbW1VQcao3SoHNz3CKD6fD5/P501rOoMQ3emIRlYFCup00i3RuaeBjgFaf5JYObhRSOC+x2pJTznEs7GMADYun+4q3fTvX+z8F0yaxu+HqqZ/on0+Wv8yKWWFrr2NaJ+PUCh1xe9GsbpVDul6hPmiZrCGmsGanBVDbGTV3oHeknQ6ma5H5UYhgffDmSOZbOLZ7JiveDaWyqNYCyGlrgup68Lvn5Q0TSZzO25oa7vJbULnn+nZX6yINHU3mQhSVRCJRMx2GeJmiM/v97N48WJXczbhN8KEl4cJTQmlTOf1HFC54FrZqOoiEYnGs4E8xLOxVB4ltnRlGF7P7QwMbMPAwJi0vSS3zJ0bJBjEs/y8orOtExYCu0HV2io6R3caz4k54DagXDECz/n9/rRKIRQK0fajNmNw8Ghyg4ORPAeUVtmIyP9LcuhYETlWVX/qsUwWS0Hwcm4nFIJNm/YEvInmGXWTEolAQ0Pu+XlJIBCAK4EOqKuvK8mhJTdKyUvF5dbgIBgMIo5LgpE2B+RmUWej8zsIOI8tvtHOBfbNn2gWS36Jzu3I2LU5V+amB5d+ZXgm+blZXFsMMln8OX369BERBjsQCFA1qgqq0i9Mra2tpbq6esTNAaXt2ajqNQAi8gTGN1qXs301cH9epXOBiBwD/BxjW/JLVb2xyCJZSoAxXznO+W9jynRu5nbcMLyXlHs0z0DA9Gj6+lKbSBeLYpg0lzJ+v59pl04jvDzMoktSWSP6aW1ttXM2aYj3jdYHTPJUmgxxrONuBT4HdAAvisjDqjqi3OdYio/bNUCZ5Nfa6m5dlJeMhF5IvvDt6cO3Z3prRDdzQJVIJsrm18ALIvIHjPnzl4B78iKVew4G3lDVVQAi8jvgeEaYrzZLaXDQNbMB8PuDnuRXCqu+k2GVUv6pNKu1TKzRbhCRx4DPOLvOVNWl+RHLNbsC62K2O4BDCnJlN5HDLJYSYnrb9GKLYHFJKBSipaWFvr4+amtrK8JqzY01mjhO2FDVfwD/SJWmwEiCfVvJISKzMcHeaG5uzrdMOfPQQw8RDAaZN29esUWpeMZ0FzC0aZEptRAKluQEg0E2bzbhQyrFas2NNdrTInKhiAyrpUWkVkQ+KyL3AKfnR7y0dAATY7YnAG/HJ1LVO1T1IFU9aPz48QUTLluWLl3KAQccUGwxLBZLkQgEAjQ0NFSU1ZqbYbRjgLMwsWsmY8x7GjCK6gngZ6rali8B0/AisJcj13rgFOCrBblyOGx+Hq64W7lyJeeffz7PP/88TU1NbNy4kYsvvtiTvC0WL8mHI1rLFtxarZXTvI4b0+cejMfn20RkFDAO2KyqG/MsW1pUdUBELgAex5g+L1DVf+X9wtFIZJEInqzgw3h8Pvnkk/n1r3/N8ccfz5IlS9h3330599xzqa+v90jwkYOdn7CUO+ms1jLxRlAKSikTazRUtR94J0+yZIWqPgo8WtCLug1KkgFPPvkk+++/P7vssgtjx45lp512or6+nsHBwZzFtZQWwekXR/8rohSWcsetR+pScZGTSVhoS5RoUBLwbMVdW1sbU6dO5ZVXXmHatGls2LCBxsZGPv74Y84880w6Ojo466yzuPPOO3nsscdQVc4666yhSUSLxTKycDuvk8hFTjHIqGdjccgkKIlLxo4dS3t7OzU1NUybNo2rrrqK888/n/Hjx9Pc3Mwll1zCXXfdxXvvvceCBQtYv349s2bNoqGhwYMbslgs5YbbeZ2oi5yoGXWxjA0yUjZOsLRqVbXNaY+Dknz961/nS1/6Eg899BDbbbcdp5xyChdeeCHd3d2sWrWKmpoaxowZw5gxY1i6dCnhcJhzzjnHk2tXItbM1zIScOONoFRc5LhWNiLybYyv1x4R+Qi4VVWzDxtpGcb222/PM888w/77789TTz1FU1MTAwMDXHTRRVx//fXcd999Qy9LTU0NV155ZbFFrgim2/kTywigFFzkuFnUeRNmIee3gX1UdYMTpfMaEblOVa/Is4ylSR5azr29vXR1ddHUZIJR1dTUsGDBAgAuvfRSwuEwF1xwAaeffjo77LCD59e3WCyWfOGmZ/MMMANj8rzE6dW0A68C54rI/5aCGXQlUFdXx6pVq5Ie9/l83HKL7UxaLJbyw806mz8AfxCRQ4H/xpg+7w9MA7YHgiIyRlX3zKukFgsQ7g0T7gkTWhcqubj3lvxjF5OWL5mYPp8P/Ab4X0xPZz/gVVWdjg2iZikAoXUh2t9rZ/XG1bQsbCG0LlRskSwWi0sy8fr8uogcgokdMx0zlPZd51hfilMtFk8IrgkSUWcR22AfwTVB27upELzusdgeUOmRqQeBPuDPzs9iKSiBSQGqpIqIRqitriUwKVBskSxljFVIhcUu6rSUDf6JfqbtOI1wT5hFX15kezUWSxlhlY2lrPDV+fDV+ayisVjKDNcGAmL4uohc6Ww3i8jB+RPNYrFYLJVCJtZotwF+4FRnuwu41XOJyoRAwBP/mxaLxTIiyETZHKKq5wM9AKr6IVCbF6lGOA899BAXXXRRscWwWCwWz8hE2fSLSDWgAI7LmkhepBrh2LDQFoul0sjEQGAe8AdgBxG5AfgK8MO8SFUG5CEqtA0LbRlOdJzWerC2VACuezaqugiziHMuxmXNCap6f74EK2WiUaFXrzZRoUMeLGSPhoX+6U9/yvjx43n++ee59tpr6enpyT1zi8ViKTKZLupcDizPkyxlQx6iQtuw0BaLpaLJxPT5JBFpdP7/oYg8JCKeTCw4ef9LRCIiclDcsctE5A0RWSEiR8fsP1BEXnWOzZNo3NMCkIeo0EnDQt93333DwkCvW7duWJjo/v7+3C9usVgseSYTA4ErVLVLRGYCRwP3ALd7JMc/gS8Dz8buFJF9gVOATwDHALc5Rgo4154N7OX8jvFIlrREo0JPngytrd7M2YwdO5bly5fT3t4+LCz0YYcdxuLFi7nrrruYNWsWEydOHAoTPW/ePEaNGpX7xS0Wh3BvmLXhtdbJqcVzMhlGi47nfAG4XVX/JCJXeyGEqr4GkKBzcjzwO1XtBVaLyBvAwSKyBhirqiHnvIXACcBjXsjjBo+jQicNCy0iw8JAx4eJtuTGTW1txRYhOfmwQklB1Kt2RCO0LGyh9bRW66nB4hmZ9GzWi8h84GTgURGpy/D8bNgVWBez3eHs29X5P35/QkRktoi8JCIvvf/++3kRNFeiYaEnTpzICy+8wA033DCkfKNhoGPDRE+fPp2gtVKqXPJhhZJmJXIir9oWi1dk0rM5GTNU9RNV3SgiOwOXuj1ZRP4K7JTg0OWq+qdkpyXYpyn2J0RV7wDuADjooIOSpsuEfNTz8WGhw+Ewl19++bAw0LFhoi0VTD6sUNJgvWrnhvUinZpM4tl8LCJPA3uJyGHObtd2uap6ZKbCYXosE2O2JwBvO/snJNhf1sSHhbZhoEcwUSuUSMQ7K5Q0w3LWq7Yln2RijXYOZgL/ceAa5+/V+RFriIeBU0SkTkQmYwwBXlDVd4AuETnUsUI7DUjWO7JYyg+vrVBcDsv9fs7rLL7+bfwdCQ9nzL5vhPnqI2vTDwOGQjB3bmHTuc3L4g2q6uoHvArUA23O9n8Av3d7fpq8v4TprfQC7wGPxxy7HHgTWAEcG7P/IIwV25vALYC4udaBBx6o8SxbtmyrfeVOJd6Tqqrv2z71fdvnWX5LfT5d6vMmv6VTfLp0iou8Dj/c/LxK54Y5c1TB/KqrzXY8S5booKARUG1oUF2yJHl+S5aYPNKkcZXfkiWqVVVGNq/SNTSY+0yWzk2aTO41H+lKDOAlzaGez2TOpkdVe0QEEalT1eUiMiVXZQegqn/AuMJJdOwG4IYE+18C9vPi+hZLxeNmWC4YRNSZEE01TxQKwcyZJq+GhuQ9L7f5uZ2fCgYharGaLl1fHwwOJk8XDMLmzenzCoVMT7Cvz5Rbsnv1Ol00bTBonpUXPVuv8sqSTJRNh4hsC/wReFJEPqQC5kmiqGoi0+uyxDRCLJYYosNy4TAsWpS4wgkEUAEUJNU8kVvlEAjQW1tFzUCEmlT5BQJGaUUr4FTpamu9Sef2mm4UV6bpvFZy6ZRIJgouj2RiIPAl59+rHUMBHwVc15JP6uvr6ezspKmpKanCWbp2KQAzmmcUUrSMUVU6Ozupr6/f6ljA+aCsyXR+2GbzAGM+HijYupiMSbc4zO/ngh9MZ/ryMLMvSaKQwL3xgt9PfXBx+srQ7zcVYCHTuc3LSwUXTeeVknOrRNwqwjzjWtk462pOBCbFnDcduNZzqQrMhAkT6OjoINUanA3/3gDAa5teK5RYWVNfX8+ECRPSJ7Qwffp0bzIKhdhz3SZEMRWAF63HIjQKlu3pY9mePmankt3vh8UulEg0rZtyKEY6t2kKrQjBnfJyq0TcKsI8k8kw2p+AMPAyZiK/Yhg1ahSTJ09OmcZ/sXmIG2/aWACJDLYnUka4nZ8oJl6+R24r/UqgWIownVJyq0TcKrg8k4mymaCqBfM/VulYRZJn3MaC8coljNv5DovFLemUUiZKpAQaB5m4m1kiIlPzJollK8LhMGvXriXk0ToAr/NzRRoXKUXFS5cwfj9vTBzNO+PqijYBaxmB+P1w2WVl8b6lVTaOG/92YCbwD8fVf3vM/hHBwDsD9LzY41lFvW84zFfXJl/sFgqFaG9vZ/Xq1bS0tOR8Xa/zqwgSWVUlIRAIDPVGk/FxQw0bmurTf/jhMKR49hZLJeJmGO24vEtR4oRCIaY+u4lAPXzvpMP40f3P4s+lJREKMa+tDQEGjziC6qef3qqCCgaDHByJEACe6+0lGAzmdM1gMEjEqVj7+vpS5+elTX6BPRdnRAYuYTzzDh3tTUUi3hkSWCxlQNqejaq+papvYVb2nwj8DPgpJv7Me/kVrzR47u4f0/oBXPcW/OW9AZ67+8c55ffWwoUAVAOR3t6h7ViOa2qiFbgOeCIS4TjHOWdCXLjdCAQCVDkR32pra5O30kMhembOZPAHP0g/tDRjBuy+e/I0mQxTpcvL4ZlFXbT96qP0vQI3vYd8BCZKRwa9KYulkshkzmYhJojZzRj3MPsAv86HUKXG3hvfpnYQahRGDZrtXHgG6AP6nd8zCdJM7eykHtP1bKiqYmpnZ+LMQiEin/oUmkY5+P1+vr7HHvx4u+34+003Je/VBIPURiJUA9rbm7wydKNI3FasbpVSKMR+nRF2+0jTpnOt5Hw+aG72RNFc/P3pXPz96akT5SPMq8VSBmSibKao6tmq+rTzmw3snS/BSok9Tj2bvmroF+ivNtu5sNdpp/E5Ea4EPl9by16nnbZ1okCAqoYGqK5G6upSLv4S4kxuExEK8avXX+c7H37I1IsvTloBv9rURA9GCW6ORHg1WY/KjSJxW7G6VUrBIFXx5sW55OcxrqJcFqM3ZbGUAJkom6Uicmh0Q0QOAf7mvUilx9QTZnPcV+q59jPVvHn/fKaeMDun/Px+Px/vvz+/nzyZucnmTqJmjdddl7pSCgToFWEAGKypyVkpPdLZSQtwJXBUVRWPJOtRRVdCV1cnVyRuK1a3SikQICIQcZGu0L2HaJTL1RtX07KwJbXC8bA3ZbGUC5kom0Mw5s9rnLDMIeDwkWKV9o8d67h5xpicFU0Un89Hc3Nzykn/wIrLCOz8eMpKKQS0AFcALaokreICASSdcsDM7bxQVcWNwD/q6pLP7bhVhm4q1uiK9DlzUufl93PMSaO57lNpzIvz0HvYZmCAHXp6kvYIbZRLiyU1mSzqtAs6C0y4N0y4J0xoXShpIKtgMMgSVZYA1YODya3MXC4A8/v9TJs2jXA4zKJFi1JbwLlZKOZ2CMvlorMQ8BxwFJAydTo/YJkQCrHnpk2mV5jEgsxGubRYUpOJI8638inISCOd54DosExEI7QsbKH1tNaECidqZRaJRFJbmYHrCt3n8+Hz+XIz784DoVCITQ9tgkFoaWmhtbW1MDI6z0owRhOSwBWNjXJpsaQmk57NiMYzh40uSTQsk6gCy6gn4vbaJWqOGwwGYRBQF2uFPLyHV5ua2AMYBfRHIrzZ1EQiVxq+Oh++Op9VNBZLAmSkxT4RkTDweswuH8bBqJvtccAHHooTf60t1DKaJvbGaVDTyUr6qEmaPn2eifa72Re7nc+ySCZPLKMhpkxMPKVUa71S5ZdJeTSMhl0bgS5gE6wH3qVY70b26TN5NxLtL81vJbv0Xr0bxfhWivVu7Kaq4zO47nByCfNZjj/gjmy3yTEsajpZvEifLE2i/W72xd1/3soiH+WR6ng25TGS3o1091vu5eH1u5HqXSn1ssjm3cjml4k1WqXwfzlu51MWL9InS5Nov5t9/5fimNd4XR6pjmdTHiPp3Ui0v5K+Fa/fjfjtciqLVGk8u48RN4yWCyLykqoeVGw5SgFbFsOx5TEcWx5bsGVhGIk9m1y4o9gClBC2LIZjy2M4tjy2YMsC27OxWCwWSwGwPRuLxWKx5B2rbCwWi8WSd6yysVgsFkvescrGI0RkHxH5hYg8ICLnFVueYiIiJ4jInSLyJxE5qtjyFBsR2V1E7hKRB4otSzEQkdEico/zTnyt2PIUm5H6PlhlA4jIAhHZICL/jNt/jIisEJE3ROT7qfJQ1ddU9VzgZKBszRw9Kos/quo3gDOAWXkUN+94VB6rVDW3IEglRobl8mXgAeed+M+CC1sAMimPSnwf3GCVjeFu4rxai0g1cCtwLLAvcKqI7CsiU0XkkbjfDs45/wksBloLK76n3I0HZeHwQ+e8cuZuvCuPSuJuXJYLMAFY5yQbLKCMheRu3JfHiMQ64gRU9VkRmRS3+2DgDVVdBSAivwOOV9W5wHFJ8nkYeFhE/gzcm0eR84YXZSEiAtwIPKaq/8izyHnFq3ej0sikXIAOjMJpo0IbuBmWx7ICi1cSVOSD94hd2dIaA/PB7JossYgERGSeiMwHHs23cAUmo7IALgSOBL4iIufmU7Aikem70SQivwBmiMhl+RauiCQrl4eAE0XkdvLv8qiUSFgeI+h9GIbt2SRHEuxLugJWVYNAMF/CFJlMy2IeMC9/4hSdTMujE6hEpRtPwnJR1U3AmYUWpgRIVh4j5X0Yhu3ZJKcDmBizPQHj0n4kYstiOLY8EmPLZTi2PGKwyiY5LwJ7ichkEakFTgEeLrJMxcKWxXBseSTGlstwbHnEYJUNICK/xYS3nyIiHSJytqoOABcAjwOvAfep6r+KKWchsGUxHFseibHlMhxbHumxjjgtFovFkndsz8ZisVgseccqG4vFYrHkHatsLBaLxZJ3rLKxWCwWS96xysZisVgseccqG4vFYrHkHatsLJYCICIXichrIrKo2LJYLMXArrOxWAqAiCwHjlXV1TH7apyFfxZLxWN7NhZLnnE8/O6OCT8RFpE7ROQJYKGIjBeRB0XkRef3aeecJhF5QkSWish8EXlLRMYV9UYslhywPRuLpQCIyBpMBNcLgC8CM1V1s4jcC9ymqotFpBl4XFX3EZF5wAeqeq2IfAF4BBivqh8U6x4sllywIQYslsLzsKpudv4/EtjXxJsDYKyINAKHYcIpo6p/FpEPCy+mxeIdVtlYLIVnU8z/VYA/RvkA4CgfO+xgqRjsnI3FUlyewAytASAi051/nwW+5uw7Ftiu4JJZLB5ilY3FUlwuAg4SkXYRWcaWCI7XAIeJyD+Ao4C1xRLQYvECayBgsZQBUQMDayBgKVdsz8ZisVgsecf2bCwWi8WSd2zPxmKxWCx5xyobi8ViseQdq2wsFovFknessrFYLBZL3rHKxmKxWCx5xyobi8ViseSd/w8m8ZCN0I1FnQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "res = ImpBI\n",
    "rho = 1e12 * np.abs(res.impedance)**2 / freq[:, None, None]\n",
    "rho_err = 1e12 * np.abs(res.error)**2 / freq[:, None, None]\n",
    "phi = np.angle(res.impedance, deg=True)\n",
    "rad_err = np.arcsin(res.error/abs(res.impedance))\n",
    "rad_err[np.isnan(rad_err)] = np.pi\n",
    "phi_err = np.rad2deg(rad_err)\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = plt.subplot(2, 1, 1)\n",
    "ax.set_xscale(\"log\", nonpositive='clip')\n",
    "ax.set_yscale(\"log\", nonpositive='clip')\n",
    "ax.errorbar(freq, rho[:,0,0], yerr=rho_err[:,0,0], fmt='k.', label=r'$\\rho_{xx}$')\n",
    "ax.errorbar(freq, rho[:,1,1], yerr=rho_err[:,1,1], fmt='g.', label=r'$\\rho_{yy}$')\n",
    "ax.errorbar(freq, rho[:,0,1], yerr=rho_err[:,0,1], fmt='r.', label=r'$\\rho_{xy}$')\n",
    "ax.errorbar(freq, rho[:,1,0], yerr=rho_err[:,1,0], fmt='b.', label=r'$\\rho_{yx}$')\n",
    "plt.xlabel('freq')\n",
    "plt.ylabel(r'apparent resistivity  $\\rho$ ($\\Omega.m$)');\n",
    "plt.legend()\n",
    "\n",
    "plt.title('Site 002 results in 2-stage RR BOUNDED INFLUENCE \\n  configuration with sites 99 and 100 as remote references')\n",
    "ax = plt.subplot(2, 1, 2)\n",
    "ax.set_xscale(\"log\", nonpositive='clip')\n",
    "ax.errorbar(freq, phi[:,0,0], yerr=phi_err[:,0,0], fmt='k.', label=r'$\\phi_{xx}$')\n",
    "ax.errorbar(freq, phi[:,1,1], yerr=phi_err[:,1,1], fmt='g.', label=r'$\\phi_{yy}$')\n",
    "ax.errorbar(freq, phi[:,0,1], yerr=phi_err[:,0,1], fmt='r.', label=r'$\\phi_{xy}$')\n",
    "ax.errorbar(freq, phi[:,1,0], yerr=phi_err[:,1,0], fmt='b.', label=r'$\\phi_{yx}$')\n",
    "plt.xlabel('freq')\n",
    "plt.ylabel(r'phase $\\phi$ (degrees)');\n",
    "plt.legend()\n",
    "plt.ylim(-180, 180)"
   ]
  }
 ],
 "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.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}