process_glider_azfp.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from echopype.convert import Convert\n",
    "from echopype import Process\n",
    "import glob\n",
    "import random\n",
    "import matplotlib.pyplot as plt\n",
    "import xarray as xr\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "from os import path\n",
    "import os\n",
    "\n",
    "\n",
    "from echopy.processing import mask_seabed as  maskSB\n",
    "from echopy.processing import mask_impulse as mIN\n",
    "from echopy.processing import mask_shoals as mSH\n",
    "\n",
    "\n",
    "\n",
    "from glider_azfp import load_glider, merge_glider_AZFP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 12min 20s, sys: 5.44 s, total: 12min 25s\n",
      "Wall time: 13min 22s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# load in the glider data\n",
    "# get from erddap first time, then save locally\n",
    "\n",
    "gdata_local = './glider_data/'\n",
    "\n",
    "# gdata_id = 'ru32-20190102T1317-profile-sci-rt'\n",
    "gdata_id = 'ru32-20190102T1317-profile-sci-delayed'\n",
    "\n",
    "# check if the glider data is already saved\n",
    "if path.exists('./glider_data/'+ gdata_id + '.nc'):\n",
    "    \n",
    "    gl = xr.open_dataset(gdata_local + gdata_id + '.nc')\n",
    "    \n",
    "# if not, load from erddap and save it    \n",
    "else:\n",
    "\n",
    "    gl = load_glider(dataset_id = gdata_id)\n",
    "    \n",
    "    if not path.isdir(gdata_local):\n",
    "        os.mkdir(gdata_local)\n",
    "    \n",
    "    # this is taking a really long time to write out (9.5 minutes)!\n",
    "    # the file is like 250MB, why so long?\n",
    "\n",
    "    gl.load().to_netcdf(gdata_local + gdata_id + '.nc')\n",
    "    \n",
    "\n",
    "# gl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# load an AZFP file\n",
    "\n",
    "you'd want to build the loop below here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "random file used:\n",
      "19011611.01A\n"
     ]
    }
   ],
   "source": [
    "path_to_raw = './data/'\n",
    "\n",
    "# alternativly, pick a random file\n",
    "all_files = glob.glob(path_to_raw+'*.01A')\n",
    "# a random file from the list (got to remove path and extention)\n",
    "fname = random.choice(all_files)[7:-4]\n",
    "\n",
    "# path to the 'fake' calibration file\n",
    "calib_file = './data/NLB_test.XML'\n",
    "\n",
    "print('random file used:')\n",
    "print(fname + '.01A')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "16:11:05 converting file 19011611.01A with NLB_test.XML, time of first ping 2019-Jan-16 11:00:01\n",
      "          ... this file has already been converted to .nc, conversion not executed.\n",
      "Initialize using average temperature recorded by instrument\n",
      "16:11:06  calibrating data in ./data_netcdfs/19011611.nc\n",
      "16:11:06  Remove noise from Sv stored in memory.\n"
     ]
    }
   ],
   "source": [
    "\n",
    "# fname = '19011711'\n",
    "# fname = '19011717'\n",
    "# fname = '19011709' #decent swarm, maybe bottom reflection?\n",
    "\n",
    "\n",
    "# converts the raw AZFP data\n",
    "# we are 'tricking' echopype with a fake .XML file\n",
    "# the real files fail because there is no temperature data\n",
    "# edit the './data/' parts so the path works for where you put the files\n",
    "data_tmp = Convert(path_to_raw + fname + '.01A', calib_file)\n",
    "\n",
    "\n",
    "# place to save processed files\n",
    "savedir = './data_netcdfs/'\n",
    "\n",
    "# converts to a netcdf we can use with xarray\n",
    "data_tmp.raw2nc(save_path = savedir )\n",
    "\n",
    "\n",
    "# read the converted data \n",
    "\n",
    "nc_path = savedir + fname + '.nc'\n",
    "\n",
    "# process so you can get calibrated, clean, data \n",
    "\n",
    "ed = Process(nc_path)   # create a processing object\n",
    "\n",
    "# need to figure out how much it matter that these vary a lot\n",
    "\n",
    "ed.temperature = gl.potential_temperature.mean().values  \n",
    "ed.salinity = gl.salinity.mean().values     \n",
    "ed.pressure = gl.pressure.mean().values     # pressure in dbar\n",
    "ed.recalculate_environment()  # recalculate related parameters\n",
    "\n",
    "ed.calibrate()           # Sv\n",
    "\n",
    "# lets remove the first few bins with ringdown?\n",
    "ed.Sv.Sv[:,:,:100] = np.NaN\n",
    "\n",
    "ed.remove_noise(SNR=2)        # denoise, this is adaptive! will change with the file! NOT WORKING NOW!\n",
    "\n",
    "\n",
    "\n",
    "azfp = ed.Sv_clean\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# merge the glider and azfp data\n",
    "\n",
    "    # count how many observation are there when the glide is pitched in the right range\n",
    "inds = (ds.m_pitch*180/np.pi < -15) & (ds.m_pitch*180/np.pi > -30)\n",
    "\n",
    "if ( inds.sum() < 3):   \n",
    "    print('')## check if glider data empty, used salinity but could use anything\n",
    "    print(' no data, skipping AZFP file')\n",
    "    print('')\n",
    "    continue \n",
    "else: # if first thing is not true, run the next line\n",
    "    ds = merge_glider_AZFP( gl, azfp, )  \n",
    "\n",
    "# ds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# processing with echopy\n",
    "\n",
    "add some processing steps from the `echopy` package\n",
    "\n",
    "try some:\n",
    "1. seabed masking\n",
    "1. impluse noise removal\n",
    "1. swarm detection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# create a mask for the bottom\n",
    "sf_mask = np.empty_like(ds.Sv)\n",
    "\n",
    "j = 0\n",
    "for f in ds.frequency:\n",
    "    sf_mask[j,:,:] = maskSB.ariza(ds.Sv.sel(frequency=f).data.T, ds.range[0,:,0].data,\n",
    "                      dc=3).T\n",
    "#     sf_mask[j,:,:] = maskSB.maxSv(ds.Sv.sel(frequency=f).data.T, ds.range[0,:,0].data.T,\n",
    "#                                   r0=0, r1=75, roff=0, thr=(-60, -80)).T\n",
    "    j+=1\n",
    "\n",
    "# create a bottom mask if any of the frequencies detect the bottom\n",
    "seafloor_mask = np.any(sf_mask, axis=0)\n",
    "\n",
    "ds.update({'seafloor_mask': (('ping_time', 'range_bin'), seafloor_mask)});"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a swarm mask\n",
    "# this takes a long time to run ...\n",
    "\n",
    "swarm_mask = np.empty_like(ds.Sv)\n",
    "\n",
    "j = 0\n",
    "for f in ds.frequency:\n",
    "    swarm_mask[j,:,:] = mSH.weill(ds.Sv.sel(frequency=f).where(~ds.seafloor_mask).data.T, \n",
    "                                  thr=-75, maxvgap=50, minvlen=0, minhlen=0)[0].T\n",
    "#     swarm_mask[j,:,:] = mSH.echoview(ds.Sv.sel(frequency=f).where(~ds.seafloor_mask).data.T,\n",
    "#                                     np.arange( len(ds.Sv.range_bin )+1),\n",
    "#                                     np.arange(len(ds.Sv.ping_time)+1),\n",
    "#                                     thr=-85,\n",
    "#                                     mincan=(5, 5),\n",
    "#                                     maxlink=(30, 2),\n",
    "#                                     minsho=(30, 2) )[0].T\n",
    "    \n",
    "\n",
    "    j+=1\n",
    "                                     \n",
    "ds.update({'swarm_mask': (('frequency','ping_time', 'range_bin'),swarm_mask)}) ;                                                                  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# plot the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Sv in swarms')"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnUAAAJcCAYAAACfVFKqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABpI0lEQVR4nO3debycdX3+/9eVBJCdIAFDEk3AAAI/RUWwblWWgtQSXLDBpVgtqS224FJl+dbWb02bWrWuSCMqdIGYsEi+iEgCAlIgGCAsYTExYQlBiEBI2LKcef/+uO9jJocz59xzZrmXuZ4+7seZuWfmvj8nwU+u+ayKCMzMzMys3EblXQAzMzMza51DnZmZmVkFONSZmZmZVYBDnZmZmVkFONSZmZmZVYBDnZmZmVkFONSZmZkBks6X9OW8y2E2Ug51ZmbWdZLeJukmSc9IekrS/0p6U97lGilJ10n6iwHnQtKr8yqT9Z4xeRfAzMx6i6RdgCuAvwLmAtsCbwc25Fkus7JzS52ZmXXbfgARcVFE9EXECxFxdUTcJWk7SWslHdz/ZknjJL0gac+BF5L0MUk3SvqqpKclrZT07rrX95Y0P20NXC7plGHKtoekBZLWS7pe0qvqrvUWSb9KWxd/Jekt6fmZJKH0O5KelfQdSTekH7szPfen6XtPScvxVFquveuuH5L+WtKy9P7/JGlfSTdLWidprqRtm//jtl7hUGdmZt32a6BP0gWS3i1pbP8LEbEBuBQ4qe79HwSuj4gnGlzvcOABYA/gK8APJCl97SJgFbA38AHgnyUdOUTZPgz8U3qtJcD/AEjaHfgp8C3g5cDXgZ9KenlEnA38EvhUROwUEZ+KiHek13tdeu7Hko4A/iX9fcYDDwFzBtz/WOCNwJuBzwOz0zJNAg4e8OdithWHOjMz66qIWAe8DQjg+8CatNVqr/QtF7J1ePlQeq6RhyLi+xHRB1xAEpj2kjQpvc8XIuLFiFgCnAd8dIhr/TQibkjD5dnAH6TX+WNgWUT8V0RsjoiLgPuBP2niV/8w8MOIuD29/pnp9SfXvedfI2JdRCwF7gGujogVEfEM8DPg9U3cz3qMQ52ZmXVdRNwXER+LiIkkLVB7A99IX74W2F7S4Wn35yHAZUNc7rd1130+fbhTes2nImJ93XsfAiYMca1H6q71LPBUep2908/WG+5aA211jfT6Tw64xuN1j18Y5PlOTdzPeoxDnZmZ5Soi7gfOJwl3RESNZALFSSStdFcMCGZZrQZ2l7Rz3blXAo8O8ZlJ/Q8k7QTsnl5nNfCqAe+tv1ZkLE/9GL0dSbpyhyqPWWYOdWZm1lWSDpD0WUkT0+eTSALcLXVvuxD4U5Iuy6G6XhuKiEeAm4B/kfQySa8FPkE6Tq6B49LlVrYlGVu3KL3OlcB+kj4kaUw68eFAklm8kLSo7TPgWgPPXQj8uaRDJG0H/HN6/QdH8vuZDeRQZ2Zm3baeZHLDIknPkYS5e4DP9r8hIhYBz5F0Wf6shXudBEwmaSW7DPiHiFgwxPsvBP6BpNv1jSShkoh4EnhPWsYnSSYxvCcifpd+7pvAB9IZuN9Kz/0jcEE6m/eDEXEN8PfAJcBjwL7A9BZ+N7OtKCJLi7GZmZmZFZlb6szMzMwqwKHOzMzMrAXpQtF3SVoi6eoBi0qfmS44/YCkYzpaDne/mpmZmY2cpF3S9ReR9LfAgRHxSUkHkiyAfRjJ+NCFwH7pmopt571fzczMrFLSHUVeS3M5Z8lIw1Z/oEvtyJYlbqYBc9LFpldKWk4S8G4eyX2GU/pQt8cee8TkyZPzLoZZ4d12222/i4hxrVxD0sG77Trq7rXP1PaMiDXtKlsvcx1mlk2TddjrXzN129veetjLMr357vs2suj2F2dLemPd6dkRMTtr+dI9gP8MeAZ4V3p6Alsv1bOK5hasbkrpQ93kyZNZvHhx3sUwKzxJA1fDb9qJx+9099hdR7HrLqOfADTsB2xYrsPMsmmyDtvmD9+yPd+dtWemN3/93KdZdPuL10fEXw5x/4XAKwZ56eyIuDzdA/hsSWcCnyJZGmewerJj495KH+rMrDskHXzUO3bgv7/7Ct787keQNM6tda379W0r8i6CWSUFQY1a5vcO+56IozLe+kLgpyShbhV1u5QAE0nWTOwIz341s0xOPH6nu//Pp3dnzBhx2ozd+LtTxz6Rd5mq4uhRJ+ZdBLPKqRFsir5MR19kC3+NSJpa9/R44P708XxguqTtJE0BpgK3tnSzITjUmdmwJB389Noab3/z9gCc9N6dWXj980hqaYyebeFgZ9Ze/S11Wf6XpaVuGLMk3SPpLuCPgNMAImIpyT7G9wJXAad2auYruPvVzDI48fid7j71z3f7/fP+1rqlD2z02Lo2OnrUiSyozcu7GGaVEBH0ZVy2rdbi8m4R8f4hXpsJzGzpBhm5pc7MhjSwla6fW+s6wy12Zu0RJF2wWY6qrNjrljozG9LAVrp+/a119/7arXXt5hY7s9YF0JcxrtUqEuvcUmdmDTVqpet30nt3ZsF1rbfWFWWLnSJxi51ZayJjK12VWuoc6sysof4Zr430t9Z9/lMtz4T9t4h4bUQcAlwBfBEg3WJnOnAQcCxwjqTRLd4rM0k/lPSEpHvqzu0uaYGkZenPsXWvtTWAOtiZjVwAfem4uuGOqsS60oe6h558nmnfuZELFz2cd1HMKmW4Vrp+7Wity7LFTkSsBPq32OmW80nCZL0zgGsiYipwTfq8YwHUwc5sZJIxddmPKih9qFv34ibuXPUMZ112N6fPuSPv4phVyZ//7Sm7DfumMWPExz+0C8C3JC2uO2Y0czNJMyU9AnyYtKWOZDudR+re1tEtdgaKiBuApwacngZckD6+ADih7nxHAqiDnVnzgqAv41GNdroKhLp6P1my2sHOrH222323bFXE2N1GA9wWEYfWHVvtmShpYbqO08BjGkBEnB0Rk4D/IdliB7q8xU5Ge0XEYwDpz/59iDIHUEkz+sPvJjZkuqmDnVlzIqAv41HLu1Zpk/LPfh3wF/GTJat5xS4v44zjXpNPecwqpJZxm53I8p4SbLHToswBNA28swF20e6Z/znxrFiz7Pq7X7O+twqq0VI3YNHAc29Y4TF2Zm1Qixp9GY5WF+4syhY7GT0uaTxA+rN/kkiZAqhZ5dUQmyLb0bEtHrqs9KFOmzYlDwb8o3LWZXc72Jm1qI8am+kb9uhrfZhxIbbYyWg+cHL6+GTg8rrzHQ+g7oY1y6YG9KFMR1Rkqc3Sd79u+9Q6tn10DRsnjEuCnbb8xZx12d0AfOjwV+ZVPLNSq2XcZqfVYcZF2WJnIEkXAe8E9pC0iqRLeBYwV9IngIeBE9NyLpXUH0A308EA6m5Ys+FFGtiyqMrs19KHOoBXXXglKz92PJvHjXWwM2uj/oU5h39fNUXESQ1eOrLB+7sWQB3szIYWAbXIFuoi4/uKrvTdr/2mnD+fMWueTp64K9asLWoZlwOoyhY7ZeOuWLPGkm3CsnW/VuWLaWVCHTjYmbVbsihnli12HOry4mBnNrik+3VUpiMqEoc6+lvkscXOlPPnM/rp9ckTBzuzltTItsVOVdZ4KisHO7OXCkQtMh55F7ZNOh1NzyeHLXb2Oe/SIYPdrCvvG8llzXpO1u11nOny52BntrVmul+rMvu1o6Euzy12hgp2596wwsHOLIPNEWzMcGxucZ06aw8HO7MtAtEXozIdWSdUFF0enchd22JnuGDnrlizofXaZthV4GBnlghEjVGZDrfUtV9TW+z07y+5DdsNedF9zrvUkyfMRqjWY10XVeFgZ+bZr93S9S12ppw/n20fXZM8cbAzy6wWGY+8C2pmNkAtsne/VuWLaR6hLpctdl514ZVDttjd9tDT7bqVWWVkbalzqCset9ZZr6shNjE609FXqI7Lkev0kiYXATcD+0talW6rMws4WtIy4Oj0eVf2eBxquZOTZt/sYGc2QNa9E2sV+ZZbNQ521ss8UaLNIuKkiBgfEdtExMSI+EFEPBkRR0bE1PTnU3XvnxkR+0bE/hHxs06UqdHkiY19wfu/d5O7Ys3qZF3jqSpb7FSRg531Kk+U6BHDrWPnFjuzxCZGszHDsbk3q5LScLCzXhQBfaFMR1W+mPZsTTxUsPvoebfkUCKz4unfEHv4lrq8S2rDcbCzXtPMNmFVGULSs6EOGge75zfVeNOXF+RUKrPiyD5RohoVYtU52FkvSbYJG5XpcPdrRTQKdmue3cjbZl2TU6nMiqGPjIOMK1Ih9gIHO+sVbqnrUfucdymj1j+fPKkLdqvWvsjRX7sun0KZFUDWgcbhqqRUHOysF9TwmLqete+589CGTS85v2zNcw521rN6bTX2XuJgZ9WXffarW+oqaMLFC5KWugGjvh3srFfVMq/x5KrEzIrFO0r0uO1Xr2Hc1TcnG8YNEuxO+M6N+RTMLCebGZVpNXYvaTJy+71xn9zu7dY6q7paOpFruKMqE/hdEw+w213LGHf1TYO22C1Z9Qx/9oNFOZXMrPuyftOtymrsvcjBzqqqj1FsjDGZjs0xOu/itoVD3SCSYDd4i90Ny37H6XPuyKdgZl3Wa6ux52VBbV6u93ewsyrKus5mLdxSV3m73bWM3W69O3kyINj9ZMlqZl15Xw6lMuuuWsaZY1UZZJwnBzuz9mpmSZOqfDF1qBvCuBtuZ6elv0meDAh2596wwsHOKq+WsVJ0qGsPBzuz9mlq8eE2DSGR9DlJIWmPunNnSlou6QFJx7TlRg041A1j/JU3sv2KR5MngwQ7d8ValdXoboVoDnZm7RI0syxT63WYpEnA0cDDdecOBKYDBwHHAudI6tgAPoe6DCZespBtH12TPHFXrPWQ7C11rkraycHOrHURXd8m7N+Bz8NWQ/SmAXMiYkNErASWA4e142aDcU2c0asuvLJhsDv3hhVcuOjhQT5lVm5ZV2P34sPt52Bn1prI2EpX11J3hKTFdceMrPeSdDzwaETcOeClCcAjdc9Xpec6YkynLlxFr7rwSh760HFsnDAuCXbakuzPuiyZVPGhw1+ZV/HM2m5zjGZTDF9N9FVkOYCiWVCbl2u4OnrUibmHS7ORCsi8MHo6hOTaiPiLRu+RtBB4xSAvnQ2cBfzRYB9rULSOcEtdk4ZqsTvrsru57aGncyiVWWdExoU7PVGic/IOVW6xs7IKsu8okaUOi4ijIuLggQewApgC3CnpQWAicLukV5C0zE2qu8xEYHXbf9lU6UNdHquxv+rCKxmzJg1vA4LdSbNvdrCzysi6TZgnSnSWg51Z85r5UtpK01lE3B0Re0bE5IiYTBLk3hARvwXmA9MlbSdpCjAVuLX1325wpQ91kE+FN+X8+Yx+en3ypC7YbewL3v+9mzzGziqhz0uaFIaDnVlzopm9Xzv0xTQilgJzgXuBq4BTI6KvIzejIqEO8qnw9jnv0kGDHSRdsQ52VnaRcTV2h7ruyDvYmZVJDbEpRmc6+jKOvcsibbH7Xd3zmRGxb0TsHxE/a9uNBlGZUAfFDHbuirUyy7oae7tCXd4Ld5ZBnsHOrXVWJs1tE1aNL6aVCnVQvGD30fNu6Xp5zNol64rs7ei6KMLCnWXRX8/lUd852FlZJGPqemv/6sqFOihWsHt+U42DvniVW+yslDbHqG52XeS+cGeZuMXObGhBE2ttdmyRke6qZKiDYgW75zb28f7v3eRgZ6VTa27hzsPLvnCnZedgZ0WX9DS4+7Uy8gp2o9Y/nzxxV6yVXNYtdtJQtygiDq07ZtdfS9JCSfcMckwjWbzzi4MUoasLd1pzHOysyJraJqwiyzJVOtRBPsFu33PnDRrsnt9U401fXtD18piNVC3jkgCRofu1DAt3llHeM2Id7KyoAjL1NNT1NpRe5UMd5BfsBuuKXfPsRt4265qul8dsJLIu3NnK3q9FWrizrBzszF6que7XauiJUAfF6opdtfZFBzsrhRoZt9hp4xpP9bq9cGeZOdiZbS1Z0sTdr5WVV4udNmx6yflVa1/k6K9d1/XymDWjlnGdp3Z+y81z4c56kiZJ+oWk+yQtlXRaen53SQskLUt/jq37TK5r6jnYmW3R3DZhDnWllEelN+HiBclXhgETJ5atec7Bzgptcw6rsRfIZuCzEfEa4M3Aqem6eWcA10TEVOCa9Hlh1tRzsDNLNLWkSd6FbZNK1sTD6Xalt/3qNYy7+ubkv7BBgp27Yq2omlm8s2oi4rGIuD19vB64j2Q5lWnABenbLgBOSB8XZk29BbV5oOr9nZg1I+vi6e5+rYBuB7vd7lrGuKtvGjTYuSvWiqrXFu5sRNJk4PXAImCviHgMkuAH7Jm+LfOaepJm9K/nt2bNmo6UeUHfjzty3SzcWmdFECE2x6hMR82hrvzyC3aDd8We8J0bu1oes+FEjy3cORhJOwGXAKdHxLqh3jrIuUHjbkTM7l/Pb9y4ce0o5qD667hRO+zQsXs04mBneQu892vPySfYDd4Vu2TVMw52Vii1Lu79WkSStiEJdP8TEZempx+XND59fTzwRHq+sGvq1Z5/Ppf7OthZnrz4cI8qUlfsklXP8Gc/WNTV8pg10uQ2YZUiScAPgPsi4ut1L80HTk4fnwxcXne+cGvqeeKE9Sq31PWwPILdbrfenTwZEOxuWPY7BzsrhMwVYjXH1L0V+ChwhKQl6XEcMAs4WtIy4Oj0eaHX1BtYv43eeeeu3t/BzvLQ3JIm1eBQV6fbwW7cDbez09LfJE8c7KyA+mIUm2P0sEdfBauSiLgxIhQRr42IQ9Ljyoh4MiKOjIip6c+n6j6Ty5p6WdTXb33r13f9/g521nUZ19nsP6qgozWxF+8c3vgrbxwy2J0+546ulsesXi3z7NdqVIhV565Y6yXeJqz9Srl4Z7cNFex+smQ1s668L4dSmUHQW4OMe4GDnfWKpsbUVaQO62ioK+vinQtq8wrVYnfuDSu4cNHDXS2PGWTfJqyKEyWqzMHOekHWJZncUjcC7Vy8sxsLd+Zh/JU3sv2KR5MnA4LdWZfd7WBnXddr+yb2Egc7q7qg9+qwroS6di/e2e2FO7tp4iUL2fbRNKgOEuzcFWvdlHnx4ap8ze0xDnZWZe5+7YCyL96ZR6X3qguvZMyap5Mng3TFOthZt9QyLt4ZUb3Zr73Cwc6qKkJsro3KdFRlslenZ7968c4RmnL+fEY/nS474DF2lpM+su2b2FeRrote5WBnVdTU7NeK9DZ0+ut1ZRfv7IZ9zru0YYudx9hZN2TvfnWos9Y42FnbRVKHZToq8sW007NfK7t4Z7dMOX/+kGPsHOyskzIv3Jl3Qa1lebfWgYOdtVdzEyWqwQNhmlS0MXYOdtZJkXFJE7fUVUMRgp1ZuzTX/VqNOsyhbgSKNsbOwc46pdcqRMs/2Lm1ztolmup+rQaHuhHKa4ydg511U9auC3e/VouDnVWDW+qsCQ52VnVZuy+qMsjYtnCws7ILsrfUVYVDXYs8K9aqrK+WbZ2nvlp1KkXbwsHOyizrmGC31NlW8hpjN1Swu+2hp7teJque7PsmVqNCtJdysLMyS8bVZTjyLmibONS1SdEmT3z0vFu6Xh6rniBj10VVakQblIOdlVFkHhNcnS+mDnVtVKQxds9vqnHQF69yi521pJlK0arNwc5Kp5nZrxX5YupQ12ZFCnbPbezj/d+7yWPsbMRqXqfO6jjYWZlEiL5atqPVOkzSP0p6dMDuWf2vnSlpuaQHJB3T8i82BIe6Dsgr2I1a/3zyxGPsrF16bI0nG56DnZVFM7Nf29RS9+/1u2cBSDoQmA4cBBwLnCNpdFvuNgiHug7Jo+Lb99x5DYOdx9jZSPTazDErBwc7yyZroBN0bgjJNGBORGyIiJXAcuCwTt3Moa6D8gp2jcbYHfKln3e9PFZu2b/lOtT1krxb68DBzobX3JImABwhaXHdMaPJW35K0l2SfihpbHpuAvBI3XtWpec6wqGuw4o0xm7tC5t505cXdL08Vl61INN4lJr7X3uOg50VXdL92tSSJtdGxKF1x+z660laKOmeQY5pwPeAfYFDgMeAr/V/rEHROsKhrqIajbFb8+xG3jbrmpxKZWXTa8sBWHMc7KzQmpr9OnwdFhFHRcTBgxyXR8TjEdEXETXg+2zpYl0FTKq7zERgddt/19SQoU7SumGO9ZJ+3anCVcWC2rzcumL1/IvJk7pgt2rti26xs2zc/WrDcLCz4upe/SVpfN3T9wL3pI/nA9MlbSdpCjAVuLXlGzYwXEvdbyJilyGOnYHnOlU4a92rv/tjtGHTS86veXajx9jZsLLvKNGaoiwHYCNThGBnNlA0ebToK5LulnQX8C7g0wARsRSYC9wLXAWcGhF9rd9ucMOFuvdnuEaW9xj5VXwTLl6w9eCB1NoXNjvY2ZCaGZPSBrkvB2Ajl3ewc2udDSZzS12LrXUR8dGI+P8i4rURcXxEPFb32syI2Dci9o+In7X8Sw1hyFAXESvqn0vaRdLu/cdg77Gh5VHxbb96DRMv/Bm8uDE54ckTllG3KsQhdHU5AGuNg50VSheb6Yoi00QJSX8p6XHgLuC29FjcyYJVWV7Bbuq3Lxp0jJ0nT1hDWQcaJxMlDi/7cgDWOgc7K4qse1dXaVxw1tmvnwMOiojJETElPfbpZMGqLq+K79Xf/bEnT1hmWbfYSZc0WVT25QCsPRzsrBACajVlOnpt79ffAM93siC9KNdg58kTlkn7ul/LsByAtY+DneUuSOqmrEcFZA11ZwI3SfoPSd/qPzpZsF7hyRNWZN3quijKcgDWXg52lqcuT/QqhKyh7j+Aa4Fb2DKm7rZOFarXePKEFVWvLQdg1eNg1+N6bLJE1lC3OSI+ExE/iogL+o+OlqzHePKEFVFknCjR+n2KsRyAtV/erXXgYNezmuhpqEimyxzqfiFphqTxA5c0sfYp4uQJB7seV8FvuZLGSnpt3uXoJQ52lpsK1mFDyRrqPkQ6rg4vadJRRZs84Vmxva2ZJQGKTNJ1/etsAncCP5L09bzL1Usc7KzbsvY0JPVXseuwrDKFurplTOoPL2nSIYWYPFFnzbMbHex6VFRnOYBdI2Id8D7gRxHxRuConMvUcxzsrOvcUreFpDcMd4Es77Hm5TXGbtzVN289ZSjlMXY9qjrLAYxJZ9h+ELgi78L0Mgc76x41cVTDcC11P0rHn+ze6AB+0I2C9qI8Kr/d7lrGxAuvHHRWrMfY9Z6sywGU4Evu/wV+DiyPiF9J2gdYlnOZepaDnXWNW+q2sitbL2Ey2PHSgVjWNrnOim0wxs7BrsdUoEKMiHnpzNq/Tp+viIj3D/UZSS+TdKukOyUtlfSl9PzukhZIWpb+HFv3mTMlLZf0gKRjOvtblZuDnXVcjwU6GCbUpduC7dNgTF3/4c21O6xoY+wc7HpHVfZNlDRO0lmSZqd7y/5Q0g+H+dgG4IiIeB3J9mXHSnozcAZwTURMBa5JnyPpQGA6cBBwLHCOpNEd+pUqwcHOOiqSccGZjoLXYVllnf1qOSvEGLs6DnY9ojrfci8n6XlYCPy07mgoEs+mT7dJjwCmAf3rdF4AnJA+ngbMiYgNEbESWM6W7c6sAQc766gea61zqCuRvMbYjbv6Jge7XpV18eHiV4g7RMQXImJuRFzSfwz3IUmjJS0BngAWRMQiYK/+xZHTn3umb58APFL38VXpucGuO0PSYkmL16xZ08KvVQ1FCHZWRZWZ6JWZQ13JFCLYefJED6lMpXiFpOOa/VBE9EXEIcBE4DBJBw/x9sH+EAaNuxExOyIOjYhDx40b12yxKinvYOfWugoKUBNHFWQOdZImSHqLpHf0H50smDWW66zYTZtf8pqDXYWVvNtC0npJ64DTSILdC5LW1Z3PJCLWAteRjJV7PF0ehfTnE+nbVgGT6j42EVjd+m/ROxzsrO1KXoc1K1Ook/SvwP8C/wf4u/T4XAfLZcPIa4zdxLlXN5w8cfTXrut6mazDSj4eJSJ2johd0p+jImL7uue7DPXZdHLFbunj7UkWK74fmA+cnL7tZJLxeqTnp0vaTtIUYCpwawd+rUpzsLO2CbL3NBS/tyGTrC11JwD7R8RxEfEn6XH8cB/ykgCdVbTJE8vWPMcJ37mx62WyDqpIhSjpJU3Jg50bYDzJvtd3Ab8iGVN3BTALOFrSMuDo9DkRsRSYC9wLXAWcGhF97fsteoeDnbVNib+UjkTWULeCZOZXs7wkQIcVYoxdnSWrnnGLXYWUffHh9Ivly4E9BiykPhnYe6jPRsRdEfH6dH27gyPi/6bnn4yIIyNiavrzqbrPzIyIfSNi/4j4WUd/uYpzsLOWNdPTUNRKrEnDbRP2bUnfAp4Hlkj6D0nf6j+Gu7iXBOiOogW7ZWue8xi7Kil3hfiXwGLgAOB2tiyafjnw3RzLZRk42FnLHOq2spikApwP/BNwE1sqxcVZbtCJJQG8HMBL5RXsdrv17uSJlzuppszdr3kXdHAR8c2ImAJ8bsCi6a+LiO/kXT4rPge7kqvA8JFmDLejxAURcQGwW//junNjh/ps3TXaviSAlwMYXB7BbtwNt7PbIge7qlIt41HQUFfnPEmfkXSppEsknS7pZXkXyoaXd2sdONiVVVPLmRS/Dssk65i6kwc597FmbuQlAbojr2A37udeoLiyMn3LLfw33QtIxup+G/gOcCDwX7mWyDJzsLMR66GuVxh+TN1Jkv4fMEXS/LrjOuDJ4S7uJQHyUbQxdg52JVad8Sj7R8QnIuIX6TED2C/vQll2DnbWtCZa6krQ25DJcC11NwFfIwliX6s7PkPS4jYcLwmQkyIGuzd9eUHXy2Qtqk6ouyOdeQ+ApMNJ1t60EnGws+Y0syRT4XsbMhluTN1DEXFdRPwB8ADJhti7AKsj4qVbC7z0814SIEdFC3Zrnt3IIV/6edfLZC2oRqADOBy4SdKDkh4Ebgb+UNLd6ZdOKwkHO2tKdb6YZpJ1R4lPkHSDvg/4AHCLpI93smDWHoUIdnXhbu0Lm91iVyYVWXyYpGdhCvCH6TEFOA54D/AnOZbLRsDBzjLpsUAH2SdKfB54fUR8LCJOBt4IfKFzxbJ2ynWv2Bc3Jifqgt2aZzd6jF2JlG3mWDru90OSduw/FxEPkUzCOiJ9/BwwKu2NeCivstrIOdjZcEQ5Z78OVodllTXUrQLW1z1fz9bryVnB5bWl2NRvX4SefzE5URfsPHmiJMr5LfdrwNuAeyXNk/QBSf9E8kX0zPQ92wL/nVcBrT0c7GxIzbTUFaseG6wOy7QEU9ZQ9yiwSNI/SvoH4BZgebru02dGVmbrtrwqwFd/98dow6aXnHewK74yzhyLiOsj4q+BfYDZwAeBs4DjSVroiIjVwM65FdLaxsHOhlS+QNeoDnti6E8lsoa63wA/YcuvfjnwGEml6IqxRPKqACdcvOAl4+vAs2ILr6Q7SqRLKL0f+CTwJuB3EfH76nsk3RpWXA521kgZv5jCoHXYBUN/IjEmy5si4kvpTXaMiOdGWkgrhgW1eV2vgLZfvYZxV9/Mmj96CxCgLQPr+2fFLvmHY7paJsuoYJXdcCT9mGS261UkCw1fD3xG0n8Au0k6Bfg48P38Smntlke9NtDRo04sRMC0VNDEJK7iTPYarA6LiFqWz2ad/foHku4F7kufv07SOSMsr/Wo4WbFermTAirReBRJb0p3qPkRsC/JsiWnA98AfghcDFwC7A98MSK+nU9JrVOKEKjyDpY2QBfrMEl/I+kBSUslfaXu/JmSlqevNWy9GKoOk7R7ljJk7X79BnAM6S4SEXEn8I6Mn7UCWlCbV7hZsV7upHhK1nXxH8CGiLgKeCvwL8B/As8AsyNiQUT8XUR8LiL8H1pFOdjZ73Wx/pL0LmAa8NqIOAj4anr+QGA6yTaFxwLnSBrd4DJD1mFZypE11BERA2e7eqeHCijarFgvd1IwJWmlS42uW8j8T0kqwfOB04D3Slo38MiroNZZDnYGaWCrZTxar8f+CpgVERsAIqJ/YsM0YE5EbIiIlcBy4LAG13hJHRYRl0TE3wOvzlKIrKHuEUlvAULStpI+R9oVa+WX66xYL3dSbCXqfgVGS+ofJ3wkcG1E7BwRu5D0MpwBTAAmkixv8uV8imm9wsGuAJqrw46QtLjumNHEnfYD3i5pkaTrJb0pPT+BrZeAW5WeG8xL6rC61zLNgcga6j4JnJoWZBVwSPrcKqKIy524KzZ/CqHa8EdBdpS4CLhe0uXAC8AvASS9Gtg+Is6JiPURsS4ivkcys8wqqgitdeBgl7vmvpReGxGH1h1bdXlKWijpnkGOaSShayzwZuDvgLmSGm0q2+hr8FB12DNZft1MoS4ifhcRH46IvSJiz4j4SEQ8meWzVh6/rwSVuVe+LRotd+K9YgugRIOMI2Im8FmSLte3pUuYQFLPrZD0YUmjJY2S9GE8hKTy8ho7PJCDXU6aGVOXoQ6LiKPSfewHHpeTNHhdGolbgRqwR3p+Ut1lJgKrG1x/qDrsb7L8ysP+6y3pXZIuTSvapZIulvTOLBe38llQmwfZZk63Tf9yJ54VW0BdGmjcpkHGRMQtEXFZ/dJLEfHr9NofBB5PjxOBD7VWaiuzMftO6dq9Xjz+MN5+wr917X6Wi58ARwBI2o9k15rfAfOB6ZK2kzQFmArc2ugijeqwiLg9SyGGDHWS/phkKYD/R1IBfhi4EvihpOOy3MAsi+FmxTrY5ag74+naMci4oYh4MCKmRcQeETEuIk6IiAf7X5d05hAft5IbrLVu829Wdu3+L5t/Ky+bf6tb7PLQvd6GHwL7SLoHmAOcnLbaLQXmAveSrDt3akR0rJdguJa6vwNOiIgfRcSdEbEkIn4InEAy0NgqKK8ui6FmxTrY5aPJJU0Oz3mQcSv8r23F1ddpr7h519zK4WDXPc3UX62OCo6IjenQtIMj4g0RcW3dazMjYt+I2D8iftbirYY0XKh7Rbom3VYi4i5gr84UyYqiaLNi176wmYO+eBW3PfR0LuXqSc19y12U8yDjVhRipod1Vn+d9ts/yDTmvGMc7Lqoi+OCi2C4UDfUlmDeLqwHFC3YPbexj/d/7yYuXPRwLuXqOW2sDDs9yLhFFanSbThFmDgBDnZd00OBDoYPdftKmj/I8f+AfbpRQOtdjYIdwFmX3e1g1wXtnDk2jJ/QhkHGLXBLXQ9xsOsRzQwfqUiwG24xu2lDvPbVdhbEiqu/AsyjAnr1d3/Mbz55IrWdd0iCnbb823vWZXcD8KHDX9n1cvWU7lR2PySZgHUPsJF0kDGwVFL/IOPNdG6QcTH+lbeuWVCb51DVCyoS1rIaMtRFxPXdKogVn0aPJvq6v7TXvufOY8VfvI++sTs72HVZ1uVKWl3SJCI2Ah9p8NpMYGYr15c0DjgFmExdvRcRH09//nMr17dyKkKwO3rUiYVpOaycjPUXFGb/6pZ1d5VZK7WrN83J7d77nHcpo59enzwZpCt21pXeta4jqjPI+HJgV2Ah8NO6w3pcEQJV3sGysoJkdG6Wo/h1WCYOddaUPCvAoYLduTescLDrhGoEOoAdIuILETE33SD7koi4JO9CWTE42FWTaHpZptJzqLOmZa0AR++6a9vvPVyw8+SJ9mqmUiy4K7xgug3Fwa6iqvPFNJMs24SdLOl2Sc+lx2JJf9aNwllxDVsBahR9z3RmLah9zruUMWueTp54VmxnVaf79TSSYPeCpHWS1ktal3ehrFgc7Cqmmfqr+HVYJsNtE/ZnwOkkG8zuTbKS++eB0xzsbMgKsMP7x045fz7bPromvZeDXcfUQBmOoleIEbFzRIyKiO0jYpf0+S55l8uKx8GuWirS05DZcC11fw28NyJ+ERHPRMTadOuL96evWY/LswJ81YVXusWu00r+LVfSAenPNwx25F0+KyYHu4pwS91L7FK/6XW/9Jy/5RqQbwU45fz5nhXbQZnH1OVd0MY+k/78Gsnamv1H/3OzQTnYlV9TEyXyLmybDBfqXhjha9Zjijwr9vQ5d+RQqooo+bfciJiRPjyOZAmTZ4C1JDtVeOKEDcnBrgJKXoc1a7hQ9xpJdw1y3A0c0I0CWnkUNdj9ZMlqt9iNVHUqxAuA1wDfAr6dPv7PXEtkpeBgV2LVqb8yG26bsNd0pRRWGXmu0L7PeZey8mPHs3ncWAbuPHHuDSt45ct39M4TTcq8o0Tni9Kq/SPidXXPfyHpztxKY9Yk7zwxMiWom9pquJa670fEQ42OrpTQSifvMXaeFdtm1fime4ekN/c/kXQ48L85lsdKZEFtXiEClVvsRqDHWuuGC3XjulIKq5y8Z8U62LVHluVMirykiaS7Jd0FHA7cJOlBSSuBm4F3ZLzGaEl3SLoifb67pAWSlqU/x9a990xJyyU9IOmYTvxOlh8Hu3JRZK/DqrKsyXDdr7tKel+jFyPi0jaXxyqkvwLMoxJ61YVXNuyKPeuyu3n4yec44ziPLhhW1m+wxa0Q39OGa5wG3MeWGf9nANdExCxJZ6TPvyDpQGA6cBDJup4LJe0XEX1tKIPZ77krtglZ66bi1mFNGa6lbleSSvFPBjnaUVmadcxQy514Vmx2yngU0VDDR7IMIZE0Efhj4Ly609NIJl6Q/jyh7vyciNgQESuB5cBhbftlrBCKEqbcYpdB1uVMKhLoYPiWuoci4uNdKYlVVp4tdvucdykr/uJ99I3d+SUtdj9ZshqAb0x/fdfLVRrlb6lr1TdIdtHZue7cXhHxGEBEPCZpz/T8BOCWuvetSs+9hKQZwAyAV77Sk3fKJs8JYfXcYpeBW+q2UtQv4FZCeVU+Q+0V6+VOhpZ18c4qkvQe4ImIuC3rRwY5N+ifTkTMjohDI+LQceM8dLmMPHmiBJpoqatK2Bku1J0s6a0DT0p6u6R9s97EA42tX16V4FCzYs+9YYUnTzTSYzPHBngrcLykB4E5wBGS/ht4XNJ4gPTnE+n7VwGT6j4/EVjdveJar3KwG0KP1WHDhbp/BtYPcv4Fkm6JrPoHGvfrH2g8Fbgmfc6AgcbHAudIGt3EfawE8gp2nhU7Aj1UGQ4UEWdGxMSImExSL10bER8h2Y3i5PRtJwOXp4/nA9MlbSdpCjAVuLXLxbYuK0JrHTjYDSZzK12F6rDhQt3kiLhr4MmIWAxMznIDDzS2wRQ12LkrdmvNLAnQQ2YBR0taBhydPicilgJzgXuBq4BTPfO1NzjYFZhb6rbysiFe2z7jPb5BMtC4vtrfaqAxUD/Q+JG69w060FjSDEmLJS1es2ZNxmJY0eQZ7BqNsTv3hhUOdnV6bTxKIxFxXUS8J338ZEQcGRFT059P1b1vZkTsGxH7R8TP8iuxdVtRgp0N0EOBDoYPdb+SdMrAk5I+AQw7eLhTA409yLg68hxjN9RyJ+6KTfXYt1yzVhQh2Lm1rk4zEyUqUocNF+pOB/5c0nWSvpYe1wN/QTJObjgeaGzDKuKsWI+xS2SuDCtSIZq1ysGuYHrsi+mQoS4iHo+ItwBfAh5Mjy9FxB9ExG+Hu7gHGltWRZwVe9Zld3PbQ0/nUKoC6aHK0KxdHOyKIRkTHNmOqEZFNlxLHQAR8YuI+HZ6XNuG+3qgsRXGUGPsPnreLYN8osc42Jk1zcGuGHqp6xUyhrp28EBjG06ei3k2GmP3/KYaB33xqt5tseux8Shm7eRgl7Nmul4rUod1LdSZFd0+5106aLB7bmMf7//eTT05xq7Xui7M2s3BLl+9NoPfoc4KJ88Wu33Ou5RR659PnnjyRKKHvuWadYKDXU7cUmdWHHlVhPueO2/IYNdLXbHuejVrDwe77su6d3WV6jCHOiu0PINdo3XsemryRI99yzXrJAe7HHSpDpP0Y0lL0uNBSUvqXuvanvYOdVZ4eXbFNpo8cciXfp5Lmbqt1xbuNOs0B7suaqL+anVMXUT8aUQcEhGHAJcAlwJd39Peoc5KoWjBbu0Lm3tnVqxb6czaysGuWyKpt7MebSBJwAeBi9JT0+jinvYOdVYaRZs80ROzYt1SZ9YRDnad11RPQ1KHHdG/r3x6zBjBbd8OPB4Ry9Lnmfa0bxeHOiuVPMfY6fkXkye9NCu2x8ajmHWTg10XNFd/Xdu/r3x6zK6/lKSFku4Z5JhW97aT2NJKBxn3tG8Xhzornbwqwld/98c9t9xJ1nXqWq2iijIexazbHOw6KEC1jEeGOiwijoqIgwc5LgeQNAZ4H/Djuo91dU97hzorpSLOiq1isOv2kgB5j0cxy4ODXYc0Eera1HZ2FHB/RKyqO9fVPe0d6qy0ijZ5AioY7Jrruji87ONRzPLiYNchmSdKtOVu09m667Xre9o71FmpFTXYzbryvhxK1X5Ndl0sKvt4FLM8Odi1VzMTJdqxTVhEfCwizh3kfNf2tHeos9IrYrA794YVnD7njhxK1QFtmihRhvEoZnlzsGuzLk32KgqHOquEIga7nyxZXfpg1+RyAK3KfTyKmSWqEuy6XIflzqHOKiPPYDdmzdPJk6oFu8yLdralRsx9PIpZESyozXOLXTsE3a7DcudQZ5WSV0U45fz5bPvomuRJlYJd1jF1tdZvVYTxKGa2tTIHu2Zm71dlAXWHOqucvILdqy68snLBLnOlmHdBzSqoCK11UOJg12Pj6cChziqqqMHuhO/cmEOpWtBU94WZtZuDXWvcUmdWEUUMdktWPcPbZl2TQ6lGptcGGZsVkYPdCEVALeNRkTrMoc4qrYjBbtXaF8sT7Hqs68KsqBzsRqjH6jCHOqu8PIPd9iseTZ6UNNh1e5swM2vMwa45CprYv7oalZhDnfWEvCrDiZcsZLdFdydPyhjseqzrwqzoHOya0MR4uqp8MXWos56RV2U47obbGffzm7ZMOqizau2LvOnLC3IpVyZdXNLEzLJxsGtCj82AdaiznpJXZbjbXcsYd/XgwW7Nsxs55Es/z6Vcw8s687UiNaJZSTjYZREosh1V4VBnPacwwa6uIln7wuZCBjvPfjUrLge7YQRQy3hUpA5zqLOelGewm3jhlfDixuRE0YNdj3VdmJWNg11jyZfO3mqtc6iznpVXZbj96jVM/fZF6PkXkxMFDnaZK0SnOrPcONgNoce+lDrUWU/LszJ89Xd/3DDYFWbyRDPdF2aWm1brstG7j21LOYoV7DKOCY7qzOB3qLOeV8Rgt+bZjYUIdlm7L8wsf63UZX1PPd22chQm2PXYcibgUGcGFDfY5d4Vm3mdugrVimYl5q7YATK31FWjDnOoM0vlHexGrX8+eVKgMXbeUcKsfBzsEsmOEhmPitRhDnVmdfKsDPc9d17xgl2Pfcs1swoJem5XHIc6swHyDnaFmjzRY99yzarCrXXQ1OLDFfli6lBnNoi8u2IbjbHr/l6xvTVzzKxKej7YBU30NFSjEnOoM2ugiMGu63vF9tgaTwNJelDS3ZKWSFqcnttd0gJJy9KfY+vef6ak5ZIekHRMfiU3SwxVj2n06K6VI5dg18ySTBWpxxzqzIaQe7DbsOkl59c8u5ED/s/PuO2h9i1B0EivdV008K6IOCQiDk2fnwFcExFTgWvS50g6EJgOHAQcC5wjqXv/apo10Kgei76+rpaj+8Eue/drVYaQdDzU+ZuulV2ewW7CxQsGnYzw4uYa7//eTcy68r7OFqAG9MXwR28tPjwNuCB9fAFwQt35ORGxISJWAsuBw7pfPLOX6tmuWHe/doS/6Vqp5bml2MQLfzboXrEA596wgj/7waKO3b/X9k0cRABXS7pN0oz03F4R8RhA+nPP9PwE4JG6z65Kz72EpBmSFktavGbNmg4V3WxrRQp2XQl3zYypq0gVllf3q7/pWumMpEJsx5iVofaKBbhh2e84+mvXtXyfQfXYt9xBvDUi3gC8GzhV0juGeK8GOTfoH0xEzI6IQyPi0HHjxrWjnGaZFCXYQRda7TymriPa/k3X33ItL81WiO0cs9JogWKAZWue68wEih77ljtQRKxOfz4BXEbyJfNxSeMB0p9PpG9fBUyq+/hEYHX3SmuWTa8EOzU1pq4alVg3Ql3bv+n6W67lqb5CHLPP5OYvoJH/327fc+cxZk06QWJAJdSRbcWyftOtRn24FUk7Stq5/zHwR8A9wHzg5PRtJwOXp4/nA9MlbSdpCjAVuLW7pTbLpleCXa8tnt7xUOdvulZF/RXi5hUPNv/haG1WwZTz57P9ikfTa21dIa19YXNbZ8YmW4D1zrfcAfYCbpR0J0k4+2lEXAXMAo6WtAw4On1ORCwF5gL3AlcBp0ZEd6cXmjVhQW1eYcJdR4Jd5p6G6gS7joY6f9O1KsuzMpx4yULG/fwm2JxmhkFmxl646OHWb9RjFWK9iFgREa9Lj4MiYmZ6/smIODIipqY/n6r7zMyI2Dci9o+In+VXerPyyTXYVUSnW+r8TdcqbahgN3r3sQ1fa4fd7lrG1H//74YTKM667O7Wg13UoJbhqFClaNZritJaB20OdrVAfdkOatWowzoa6vxN13pBowqx76n2dIEOp9HuE5AEu9Pn3DHyi/fYzDGzXtVfj2nMNoO+Pmbv8V0rS1uDXY/1NnhHCbM2yPub7lAzY3+yZDUnfOfGkV0463i6ilSIZr1sQW0esfmlu9gAbF79WNfKMfo1+7Un2AVJC1yWoyJVmEOdWZv0B7sx4/bI5f77njuP0U+vT54MCFlLVj3DNi+fdMCILtwj33DNLP8vqAB99/0aaEeLXTOtdK3VY5IOkXRL/+5Zkg6re61rO2U51Jm10YLaPDav+V1u99/nvEsbLnmibbbbsekLZv2WW5HxKGY2dLAb2D07euedO1qWloJd0M21Nr8CfCkiDgG+mD7v+k5ZDnVmbZb3N92XLHnSih77lmtmiUb12MDu2b716ztelpEHu+611KUX2CV9vCtblmObRhd3ynKoM+uAvIPdxEsWstPS3yRPWgl2PfYt18y2yFKPjXnlxC6UZITBrvkxdUf071aVHjOGvsFWTgf+TdIjwFeBM9PzmfeEbgeHOrMOyTvYjb/yxtaDXS2grzb8UWttQWUK8i3XzLbWqB4bveuuAGx+eBWjd9utiyVqRiTLMmU6AuDa/t2q0mN2/dUkLZR0zyDHNOCvgE9HxCTg08AP+j82eME6w6HOrIN+XyEO2BpszNR9f39+zMtfzpi99qQTxl95I+N+fhOj164nan2bm75AcxXi4WX/lmtmL7WgNo9R271sq3N9zzyz5fHatV0uUUbNjKnLkLMi4qiIOHiQ43KSjRQuTd86jy1fPru6U5ZDnVmHLajNY+DWYJuX9beg1dj85JNsfvyJQT7ZHrvdtYx9vn8pm55YeWfzn26qQlxU9m+5Zja4n7/wXy19vpvr3P1eNDHRq+XOBlYDf5g+PgJYlj7u6k5ZYzp1YTPbYkFtXmc3re6UGtlmtmaoECPiqEavSfpP4LT06TzgvPSx94M2K4hW6rHNqx9jzN7ju7reHdDE0JOWvyueAnxT0hjgRWBGcvtYKql/p6zNdHinLLfUmXVJ3mPsRqR7M8cK8S3XzIY20nps1Pbb5xPosmxz2IatDiPixoh4Y7qD1uERcVvda13bKcuhzqyLRhzslNP/Vbs3+/UU4GvpPtH/TN23XLwftFmhjKQeq73wQgdKMowuhrqicKgz67IRBbtofcDHyDSzzlMLdynIt1wzy6Y0PQ+Zv5g61JnZCC2ozStHpVjL+E23IhWimWVX/DqstwIdONSZ5arwlWJk7brIqyXRzPJU6DqsmcWHK7LVoUOdWc4KXSk2txq7mfWgBbV5vx/3O2bPcSO6Rkf2kI0gopbpqEol5lBnVgCFDXYZK0W31Jn1tgV9PwZg8xNrRvT5ju0hm/mLqUOdmbVRIYNd1sU7q1EfmlkLCleHNTNJoiJ1mEOdWYGUulI0s55XqDos8JImZpavoSpFjdmmiyXBs1/NrGnFCXbNfCmtRh3mUGdWQI0qxdi8qcslcUudmTWvCMEuIohaLdtRkTrMoc6soApRKfbViL6+4Y+aJ0qY2daKUId5TJ2ZFUbulWLmiRIVqRHNrK1yrcNqAX192Y6KfDF1qDMruAW1eYzeaad8bt6/XMmwh0OdmQ0uv2AXRC3jUZE6zKHOrASuWvejXO6bZLbeqRDNrDNyCXYR2b+YVqT/1aHOrCTy+bbrljoza4+u12GR8UtphdbadKgzK5FuV4rZK8SK1Ihm1lELavPQ6NHduNULG3kxc0vdxngR4IVuFKyTHOrMSub3wU4d/7/vXU/FbzNViE/G4wB3dbpAZlZ+V2+a043b3L2RDTwfzw77xs2xiSf5LcCVHS9VhznUmZXQgto8urDf6gVP8CibY+i18Z6PZ3mB9QALOl0gM6uGTvc6RERM4QBWct+w732Y5UxkHyJiQ0cL1QUOdWYl1elujIjYMJF9eJhlQ75vJfczhdcQni1hZk3odLC7gxtHvcBzQ7bWbY5NPMEq7ueOl3W0MF3iUGdWYldvmtPRivF+7njZUK11/a10d3Cj6xIza1qj+mv07mNbvnaW1roqtdKBQ51ZJXQq2A3XWudWOjNr1WD1V99TT7fl2kO11lWtlQ4c6swqY6hgt+yHh474uo1a69xKZ2btUl9/jdlrz9+fW1CbxwsnHM60e58c0XWHaq2rWisdONSZVcrAYPf45a/hwlU3seLYH4z4mo1a69xKZ2bt1F9/bX78ia3qshsv/VxL1x2sta6KrXTgUGdWOQtq81j9+bfws9V3sORNydIBT9eeb+maA1vr3EpnZp3Q3zo30KcOuHbE1xysta6KrXTgUGdWSUtnfZpR6f+9Xz5qR9626JSWrjewtc6tdGZWJvWtdVVtpQOHOrPKGvWKX/No33qmXvcxJn2l9f+r97fWrYun3UpnZqVS31pX1VY6cKgzq7RJEx6jb8NoHn3HTi1fq7+17g5+6VY6Myud/ta6x3mkkq10ACp7vSxpPfBA3uUYgT2A3+VdiBFwuburneV+VUSMa+UCkrYDTgG+61DXHq7DuqqMZQaXu1876rDXA2MjYuSD9AqsCqFucUSMfL2GnLjc3eVyW1GV9e+4jOUuY5nB5bbs3P1qZmZmVgEOdWZmZmYVUIVQNzvvAoyQy91dLrcVVVn/jstY7jKWGVxuy6j0Y+rMzMzMrBotdWZmZmY9z6HOzMzMrAJKHeokHSvpAUnLJZ2Rd3n6SZok6ReS7pO0VNJp6fndJS2QtCz9ObbuM2emv8cDko7Jr/QgabSkOyRdkT4vfLkl7SbpYkn3p3/uf1CScn86/W/kHkkXSXpZGcptrStq/QXlrsPKWH+lZSldHeb6q4AiopQHMBr4DbAPsC1wJ3Bg3uVKyzYeeEP6eGfg18CBwFeAM9LzZwD/mj4+MC3/dsCU9PcanWP5PwNcCFyRPi98uYELgL9IH28L7Fb0cgMTgJXA9unzucDHil5uH235uy9s/ZWWr7R1WBnrr7Q8parDXH8V8yhzS91hwPKIWBERG4E5wLScywRARDwWEbenj9cD95H8H2Aayf9xSX+ekD6eBsyJiA0RsRJYTvL7dZ2kicAfA+fVnS50uSXtArwD+AFARGyMiLUUvNypMcD2ksYAOwCrKUe5rTWFrb+gvHVYGesvKHUd5vqrYMoc6iYAj9Q9X5WeKxRJk4HXA4uAvSLiMUgqTWDP9G1F+l2+AXweqNWdK3q59wHWAD9Ku13Ok7QjBS93RDwKfBV4GHgMeCYirqbg5ba2KM3fZcnqsG9QvvoLSliHuf4qpjKHOg1yrlDrs0jaCbgEOD0i1g311kHOdf13kfQe4ImIuC3rRwY5l8ffwRjgDcD3IuL1wHMkzf6NFKLc6ViTaSRdEXsDO0r6yFAfGeRcof6bt8xK8XdZpjqsxPUXlLAOc/1VTGUOdauASXXPJ5I0/RaCpG1IKsP/iYhL09OPSxqfvj4eeCI9X5Tf5a3A8ZIeJOkOOkLSf1P8cq8CVkXEovT5xSQVZNHLfRSwMiLWRMQm4FLgLRS/3Na6wv9dlrAOK2v91V+WstVhrr8KqMyh7lfAVElTJG0LTAfm51wmACSJZGzEfRHx9bqX5gMnp49PBi6vOz9d0naSpgBTgVu7Vd5+EXFmREyMiMkkf57XRsRHKH65fws8Imn/9NSRwL0UvNwk3RZvlrRD+t/MkSRjl4pebmtdYesvKGcdVtb6C0pbh7n+KqK8Z2q0cgDHkczK+g1wdt7lqSvX20iale8ClqTHccDLgWuAZenP3es+c3b6ezwAvLsAv8M72TJ7rPDlBg4BFqd/5j8Bxpak3F8C7gfuAf6LZGZY4cvtoy1/94Wsv9KylboOK1v9lZaldHWY66/iHd4mzMzMzKwCytz9amZmZmYphzozMzOzCnCoMzMzM6sAhzozMzOzCnCoMzMzM6uAMXkXwIpH0j8CzwK7ADdExMIG7zsB+HVE3Nu90pmZNeb6y3qZW+qsoYj4YqMKMXUCcGCXimNmlpnrL+tFDnUGgKSzJT0gaSGwf3rufEkfSB/PknSvpLskfVXSW4DjgX+TtETSvpJOkfQrSXdKukTSDnXX+ZakmySt6L9m+trnJd2dfmZWem5fSVdJuk3SLyUd0PU/EDMrDddfZgl3vxqS3kiyrc7rSf6buB24re713YH3AgdEREjaLSLWSppPsmr7xen71kbE99PHXwY+AXw7vcx4klXqDyDZLuZiSe8m+bZ8eEQ8n94HYDbwyYhYJulw4BzgiM79CZhZWbn+MtvCoc4A3g5cFhHPA6SVXb11wIvAeZJ+ClzR4DoHp5XhbsBOwM/rXvtJRNSAeyXtlZ47CvhR/30j4ilJO5FsCj0v2U4QSLaeMTMbjOsvs5RDnfVruF9cRGyWdBjJhs3TgU8x+DfP84ETIuJOSR8j2X+x34a6x6r7OfC+o4C1EXFIE2U3s97m+ssMj6mzxA3AeyVtL2ln4E/qX0y/fe4aEVcCp5NsPA2wHti57q07A49J2gb4cIb7Xg18vG7syu4RsQ5YKenE9JwkvW7Ev5mZVZ3rL7OUQ50REbcDPwaWAJcAvxzwlp2BKyTdBVwPfDo9Pwf4O0l3SNoX+HtgEbAAuD/Dfa8iGZ+yWNIS4HPpSx8GPiHpTmApMG3Ev5yZVZrrL7MtFNGw1drMzMzMSsItdWZmZmYV4FBnZmZmVgEOdWZmZmYV4FBnZmZmVgEOdYakD0u6Ou9ymJmZ2ch59muJSXob8BXgIKAPuA84PSJ+lWvBzMzMrOu8o0RJSdqFZLubvwLmAtuSbJezYajPlZmSfXeUbtdjZmZmddz9Wl77AUTERRHRFxEvRMTVEXGXpO0krZV0cP+bJY2T9IKkPQdeSNLHJN1Y9zwkfVLSMklPS/qu6jYyHPDZwyQtlrRO0uOSvp6ev0DSZ9PHE9Jr/nX6/NWSnkpXWx8r6QpJa9J7XSFpYt31r5M0U9L/As8D+/RfKy3fekn/JGlfSTen5Zgradv083uk11yb3vOXkvzfvZmZVY7/cSuvXwN9aXh6t6Sx/S9ExAbgUuCkuvd/ELg+Ip7IeP33AG8CXpd+9pgG7/sm8M2I2AXYl6TVEJKV29+ZPv5DYEX6E+AdwC8j6fsfBfwIeBXwSuAF4DsD7vFRYAbJyvAPpeeOBd4IvBn4PDCbZCX3ScDBbPndPwusAsYBewFnMcQ+kWZmZmXlUFdS6R6DbyMJKN8H1kiaL2mv9C0XsnWo+1B6LqtZEbE2Ih4GfsGW/RIH2gS8WtIeEfFsRNySnr8eeHvaKvYOkrF/b01f+8P0dSLiyYi4JCKej4j1wEy2hL9+50fE0ojYHBGb0nP/GhHrImIpcA9wdUSsiIhngJ8Br68r33jgVRGxKSL6w6SZmVmlONSVWETcFxEfi4iJJK1TewPfSF++Fthe0uGSXkUSyi5r4vK/rXv8PLBTg/d9gqQr+H5Jv5L0nrRsvwGeTe/7dpLxf6sl7U9dqJO0g6T/kPSQpHUkm3PvJml03T0eGeS+j9c9fmGQ5/3l/TdgOXC1pBWSzhj61zYzMysnh7qKiIj7gfNJwh3pZIK5JK11HwKuSFvC2n3fZRFxErAn8K/AxZJ2TF++HvgAsG1EPJo+/zNgLMnm25B0j+4PHJ524b4jPV8/hm/ELWsRsT4iPhsR+wB/AnxG0pEjvZ6ZmVlROdSVlKQDJH22f1KBpEkkAe6WurddCPwpyVizZrpemynHRySNS0Pk2vR0X/rzeuBTJK1vANcBfwPcGBH979mZpGVtraTdgX9oc/nek07MELAuLVvfMB8zMzMrHYe68loPHA4skvQcSZi7h6TlC4CIWAQ8R9It+7MOleNYYKmkZ0kmTUyPiBfT164nCW39oe5GYIe655B0F28P/C79Ha5qc/mmAgtJuoJvBs6JiOvafA8zM7PcefFhMzMzswpwS52ZmZlZBTjUmVkhpItI3yVpiaSrJe1d99qZkpZLekBSozUTzcx6mrtfzawQJO2Srr+IpL8FDoyIT0o6ELgIOIxkfOhCYL+6yTZmZob3fjWzBiS9Dzi5iY88ERGnjPR+/YEutSNblrKZBsxJd0pZKWk5ScC7eaT3MjOrotKHuj322CMmT56cdzHMCu+22277XUSMa+IjR1z6o/HHv/G122V68x9/eDWSfkWypVu/2RExO+sNJc0kWcvwGeBd6ekJbL1Uz6r0nJmZ1Sl9qJs8eTKLFy/OuxhmhSfpoeHftbU99hjF3nuPHv6NwDbbQBrgGoY4SQuBVwzy0tkRcXlEnA2cLelMkjUO/4GtF6Lu53EjZmYDlD7UmVnn1KJGX9Tadr2IOCrjWy8EfkoS6lYBk+pemwisbluhzMwqwrNfzayhIKhlPFptOpM0te7p8cD96eP5wHRJ20maQrKg9K0t3s7MrHLcUmdmDSWBLWtLXcs9orMk7Q/UgIeATwJExFJJc4F7gc3AqZ75amb2Ug51ZtZQLYK+jMsetRrpIuL9Q7w2E5jZ4i3MzCrNoc7MGurvWjUzs+JzqDOzhmoEfRlDnaOfmVm+HOrMrKHNBJsyzn5tfaqEmZm1oqOzXyX9UNITku6pO7e7pAWSlqU/x9a95v0dzQqk1sRhZmb56vSSJucDxw44dwZwTURMBa5Jn5Pu7zgdOCj9zDmSsq16amYd0d/9muVwO52ZWb46Guoi4gbgqQGnpwEXpI8vAE6oOz8nIjZExEqgf39HM8tJLaAv42FmZvnKY/HhvSLiMYD0557p+QnAI3Xva7i/o6QZkhZLWrxmzZqOFtaslzXT/epcZ2aWryLtKJF5f8eImB0Rh0bEoePGNbM/uZk1o4boy3iYmVm+8gh1j0saD5D+fCI97/0dzQqmFtmPjGsUm5lZh+QR6uYDJ6ePTwYurzvv/R3NCqS5ljq31pmZ5amj69RJugh4J7CHpFXAPwCzgLmSPgE8DJwI3t/RrIhqkLlr1Q11Zmb56mioi4iTGrx0ZIP3e39HswKphaiFW+DMzMrAO0qYWUO1JiZBuKXOzCxfDnVm1lAfo+jLPPTWLXpmZnlyqDOzhvpCbIpsoc4tdWZm+XKoM7OG+hhFn0OdmVkpONSZWUOBqLn71cysFBzqzKyhZnaLcEudmVm+HOrMrKFauPvVzKwsHOrMrKEaybIm2bj71cwsTw51ZtZQrYklTdxSZ2aWL4c6M2uoFsrc/epUZ2aWL4c6M2uoxqjMs1/D3a9mZrlyqDOzhpKWOoc1M7MycKgzs4aSvV89ps7MrAwc6sysoc2MZlNkqybc/Wpmli+HOjNrqM/dr2ZmpeFQZ2YNhSdKmJmVhkOdmTXU18ySJmZmliuHOjNrKFDmHSU8UcLMLF8OdWbWUDN7v+Kxd2ZmuXKoM7OG+rykiZlZaTjUmVlDEaLmFjgzs1LwCGgza6iPUZmPds1+lfQ5SSFpj7pzZ0paLukBSce05UZmZhXjljozaygQtS7OfpU0CTgaeLju3IHAdOAgYG9goaT9IqKvawUzMysBt9SZWUPJmLpsR5vG1P078Hm2HqI3DZgTERsiYiWwHDisPbczM6sOt9SZWUN9MYpNMTrTewMhaQYwo+707IiYneXzko4HHo2IO6WtunInALfUPV+VnjMzszoOdWbWUC1GNdX9mga4hiFO0kLgFYO8dDZwFvBHg31ssFtlLpSZWY9wqDOzhmpN7CgRGWbJRsRRg52X9P8BU4D+VrqJwO2SDiNpmZtU9/aJwOpMhTIz6yEeU2dmDdXSHSWyHK2IiLsjYs+ImBwRk0mC3Bsi4rfAfGC6pO0kTQGmAre2+ruZmVWNW+rMrKEa2XeU6FR/aEQslTQXuBfYDJzqma9mZi/lUGdmDdWCJhYfbt8ixWlrXf3zmcDMtt3AzKyCHOrMrKFaurBwFp65YGaWL4c6M2uo5m3CzMxKw6HOzBoKRC1zS53Dn5lZnhzqzKyhvhB9bqkzMysFhzozayia6H4ND6ozM8uVQ52ZNVRDTe0oYWZm+XGoM7OGNscoNmVep87dtGZmeXKoM7OGosm9X83MLD8dra0lTZL0C0n3SVoq6bT0/O6SFkhalv4cW/eZMyUtl/SApGM6WT4zG1qN7FuFeUidmVm+Ov0VfDPw2Yh4DfBm4FRJBwJnANdExFTgmvQ56WvTgYOAY4FzJI3ucBnNrIFaOvs1y9HOHSXMzKx5HQ11EfFYRNyePl4P3AdMAKYBF6RvuwA4IX08DZgTERsiYiWwHDisk2U0s8aCpPs1y+HZr2Zm+eraYBlJk4HXA4uAvSLiMUiCH7Bn+rYJwCN1H1uVnht4rRmSFktavGbNmo6W26yX9e/9muUwM7N8dSXUSdoJuAQ4PSLWDfXWQc695Pt/RMyOiEMj4tBx48a1q5hmNkDW8XQeU2dmlr+OhzpJ25AEuv+JiEvT049LGp++Ph54Ij2/CphU9/GJwOpOl7FTjnzXv+RdBLOWRMZWuprH1JmZ5a7Ts18F/AC4LyK+XvfSfODk9PHJwOV156dL2k7SFGAqcGsny9hJ1/zizLyLYNaSJLBlG1NnZmb56vQ6dW8FPgrcLWlJeu4sYBYwV9IngIeBEwEiYqmkucC9JDNnT42Ivg6X0cwaCJrYJqzDZTEzs6F1NNRFxI007pM5ssFnZgIzO1YoM8usf7ycmZkVn3eUMLOG+mqj2FzLuE2YZ8CameXKoc7MGop0SRMzMys+hzoza6iZNeg8ps7MLF8OdWbWkMfUmZmVh0OdmTXU3OxXhz8zszw51JlZQ7VmxtS5/9XMLFcOdWbWUHhfVzOz0nCoM7OGPFHCzKw8HOrMrKFAXn/OzKwkHOrMrKFoavarw5+ZWZ4c6sysoWYmSrj71cwsXw51ZtZYNNH96lRnZpYrhzoza6ivNoq+rHu/drgsZmY2NIc6M2soaGbvV4+pMzPLk0OdmTUUkRxmZlZ8DnVm1lAze786+5mZ5cuhzswa80QJM7PSyDYC2sx6Uv+OElmOVkn6R0mPSlqSHsfVvXampOWSHpB0TMs3MzOrILfUmVlDQdfH1P17RHy1/oSkA4HpwEHA3sBCSftFRF9XS2ZmVnBuqTOzhiLtfs10dK4Y04A5EbEhIlYCy4HDOnc7M7Nycqgzs8Yie7ADIWmGpMV1x4wm7/gpSXdJ+qGksem5CcAjde9ZlZ4zM7M67n41s4aaGi8XEBGzgdmN3iJpIfCKQV46G/ge8E/Jlfgn4GvAxxl8ATxPyzAzG8Chzswaavc6dRFxVJb3Sfo+cEX6dBUwqe7licDq9pXKzKwa3P1qZg1FjKJWy3a0mv0kja97+l7gnvTxfGC6pO0kTQGmAre2eDszs8pxS52ZNRR0tZ/zK5IOSW/5IPCXABGxVNJc4F5gM3CqZ76amb2UQ52ZNRTpRInu3Cs+OsRrM4GZXSmImVlJOdSZWWPNNNV56oKZWa4c6sysoaCJbcIy7hFrZmad4VBnZo01MfvVDXVmZvlyqDOzhrYsLGxmZkXnUGdmjYWSw8zMCs+hzswaamrxYfe/mpnlyqHOzIbmsGZmVgoOdWbWkMfUmZmVx4hCnaSxwKSIuKvN5TGzIqlB1BzqzMzKIPPer5Kuk7SLpN2BO4EfSfp654pmZoUQTRxmZpabzKEO2DUi1gHvA34UEW8EjupMscysGNTkYWZmeWkm1I2RNB74IHBFh8pjZkXiVjozs9JoJtT9X+DnwPKI+JWkfYBlQ31A0ssk3SrpTklLJX0pPb+7pAWSlqU/x9Z95kxJyyU9IOmYkfxSZtYmDnVmZqWROdRFxLyIeG1E/HX6fEVEvH+Yj20AjoiI1wGHAMdKejNwBnBNREwFrkmfI+lAYDpwEHAscI6k0U3+TmbWLv2LD2c5zMwsV5lnv0oaB5wCTK7/XER8vNFnIiKAZ9On26RHANOAd6bnLwCuA76Qnp8TERuAlZKWA4cBN2ctp5m1TzOLD7uxzswsX80saXI58EtgIdCX9UNpS9ttwKuB70bEIkl7RcRjABHxmKQ907dPAG6p+/iq9NzAa84AZgC88pWvbOJXMLOmeUcJM7NSaCbU7RARX2j2BhHRBxwiaTfgMkkHD/H2wfpwXvJPRUTMBmYDHHroof6nxKxT3LVqZlYazUyUuELScSO9UUSsJelmPRZ4PJ1JS/rzifRtq4BJdR+bCKwe6T3NrDWKJo68C2tm1uOGDXWS1ktaB5xGEuxekLSu7vxQnx2XttAhaXuSde3uB+YDJ6dvO5mka5f0/HRJ20maAkwFbh3B72Vm7eDZr2ZmpTFs92tE7NzC9ccDF6Tj6kYBcyPiCkk3A3MlfQJ4GDgxvddSSXOBe4HNwKlp962Z5cXdr2ZmpdDM7NdrIuLI4c7VS/eGff0g558EBv1cRMwEZmYtl5l1UC09snBrnZlZroYNdZJeBuwI7JEuEtz/tX0XYO8Ols3M8uauVTOz0sjSUveXwOkkAe72uvPrgO92oExmVhSe/WpmVhpZxtR9E/impL+JiG93oUxmViByS52ZWSk0s07deZI+A7yNpEPml8C5EfFiR0pmZvlrpvvV4c/MLFfNhLoLgPVAf2vdScB/kc5cNTMzM7P8NBPq9o+I19U9/4WkO9tdIDMrjv6Fhc3MrPia2VHiDklv7n8i6XDgf9tfJDMrjP6JElkO7ylhZparZlrqDgf+TNLD6fNXAvdJuhuIiHht20tnZvlzS52ZWSk0E+qO7VgpzKyYSrhOnaT5wBzg8oh4Lu/ymJl1S+bu14h4CJgEHJE+fg4YFREPpc/NrGL6x9RlOQrkaySz9O+VNE/SB9JF1M3MKq2ZbcL+ATgU2B/4EbAt8N/AWztTNDPLXQm3CYuI64Hr0z2njwBOAX5IsguOmVllNdP9+l6SfVxvB4iI1ZJ27kipzKwQROFa4TKRtD3wJ8CfAm8gWZLJzKzSmgl1GyMipKSKl7Rjh8pkZkVRwm3CJP2YZGLXVcB3gOsjImt7o5lZaTWzpMlcSf8B7CbpFGAh8P3OFMvMCiGaONpA0t9IekDSUklfqTt/pqTl6WvHNPjsmySNJxkesi9wM8m+1d+QtHt7SmhmVlyZW+oi4quSjgbWkYyr+2JELOhYycwsf01Mgmi1m1bSu4BpwGsjYoOkPdPzBwLTgYOAvYGFkvaLiL4Bl/gP4KiIuErSO4B/Af4GOASYDXygtRKamRVbM92vpCHOQc6sl3RvTN1fAbMiYgNARDyRnp8GzEnPr5S0HDiMpCWu3uiIeCp9/KfA7Ii4BLhE0pKOl97MLGfDdr9KWi9pXaOjG4U0s3w0u6SJpBmSFtcdM5q43X7A2yUtknS9pDel5ycAj9S9b1V6bqDRkvq/qB4JXFv3WlNfYM3MymjYii4idgaQ9H+B3wL/RTIp7sOAZ7+aVVmT4+UiYjZJV+egJC0EXjHIS2eT1EdjgTcDbyIZx7sPg+8/NlipLiJZyuR3wAvAL9N7vhp4JvtvYWZWTs18ez0mIg6ve/49SYuArzT6gJmVXJt3lIiIoxq9JumvgEsjIoBbJdWAPUha5ibVvXUisHqQa8+UdA0wHrg6vQ4kPRJ/06ZfwcyssJqZ/don6cOSRksaJenDwMCBymZWIU11v7Ye/n5CslgwkvYjWeD8d8B8YLqk7SRNAaYCtw52gYi4JSIuq98eLCJ+HRG3t1w6M7OCaybUfQj4IPB4epyYnjPruHf90b/mXQTrvB8C+0i6h2Tv1pMjsRSYC9xLsvbcqYPMfDUz63nNLGnyIMkstEFJOjMi/qUdhTIb6BdXfyHvIvQkBahLy/ZGxEbgIw1emwnM7E5JzMzKqZmWuuGc2MZrWcEde8gX8y6CdUOXFx82M7ORa+c0/3LtJWQt2Thuh7yLYN3gwGZmVhrtDHWu+nvItQvOyLsI1gWiiZ0iXAOYmeXKLXVm1phb6szMSqOdoW5eG69lZkXQzN6vnS2JmZkNI3OokzQOOAWYXP+5iPh4+vOf2104M8uZW+rMzEqjmZa6y0m23VmIFx026w0OdWZmpdFMqNshIrxYmFkP8UQJM7PyaGaduiskHdexkphZ8XidOjOz0mgm1J1GEuxekLRO0npJ6zpVMDMrAIc6M7PSaGabsJ07WRAzKx41MfvVzMzyNWyok3RARNwv6Q2DvR4Rt7e/WGZWBE3t/erwZ2aWqywtdZ8BZgBfY+tqW+nzIzpQLjMrAnetmpmVxrBj6iJiRvrwOOCnwDPAWmB+es7Mqspj6szMSqOZJU0uANYB30qfnwT8J/DBdhfKzIrDO0WYmZVDM6Fu/4h4Xd3zX0i6s90FMrMCcSucmVlpNLOkyR2S3tz/RNLhwP+2v0hmVhT9iw9nOvIurJlZjxs21Em6W9JdwOHATZIelLQSuBl4R5abSBot6Q5JV6TPd5e0QNKy9OfYuveeKWm5pAckHTOyX8vM2qKZMXVu0TMzy1WW7tf3tOE+pwH3Abukz88AromIWZLOSJ9/QdKBwHTgIGBvYKGk/SLCe82a5cFhzcysNLLMfn1oqGO4z0uaCPwxcF7d6WkkEy9If55Qd35ORGyIiJXAcuCwJn4fM2ujzF2vDn5mZrlrZkzdSH0D+DxQv4TpXhHxGED6c8/0/ATgkbr3rUrPbUXSDEmLJS1es2ZNRwptZil3vZqZlUJHQ52k9wBPRMRtWT8yyLmX/HMREbMj4tCIOHTcuHEtldHMGnNLnZlZeTSzpMlIvBU4XtJxwMuAXST9N/C4pPER8Zik8cAT6ftXAZPqPj8RWN3hMppZA6p5mzAzs7LoaEtdRJwZERMjYjLJBIhrI+IjJLtRnJy+7WTg8vTxfGC6pO0kTQGmArd2soxmNgTPfDUzK41Ot9Q1MguYK+kTwMPAiQARsVTSXOBeYDNwqme+muXMgc3MrBS6Fuoi4jrguvTxk8CRDd43E5jZrXKZ2RA8Xs7MrDTyaqkzszJopmvV4c/MLFcOdWbWkCJQZEtr3ibMzCxfDnVmNjS3wJmZlYJDnZk15DXozMzKw6HOzBrzciVmZqXhUGdmDTXVUufwZ2aWK4c6Mxuaw5qZWSk41JlZYx5TZ2ZWGh3dJszMyq1/79csR8v3kn4saUl6PChpSd1rZ0paLukBSce0fjczs+pxS52ZNRaRHFnf29Kt4k/7H0v6GvBM+vhAkr2jDwL2BhZK2s9bCJqZbc0tdWbWkNgyWWLYo133lAR8ELgoPTUNmBMRGyJiJbAcOKxNtzMzqwyHOjNrLJo4AEkzJC2uO2aM4K5vBx6PiGXp8wnAI3Wvr0rPmZlZHXe/mllDzY6Xi4jZwOyG15MWAq8Y5KWzI+Ly9PFJbGmlg8EbAT19w8xsAIc6MxtaG9epi4ijhnpd0hjgfcAb606vAibVPZ8IrM5YKjOznuHuVzNrKPN4uva1mx0F3B8Rq+rOzQemS9pO0hRgKnBr2+5oZlYRbqkzs8aamf3aHtPZuuuViFgqaS5wL7AZONUzX83MXsqhzswa6/LiwxHxsQbnZwIzu1cSM7Pycagzs8bqZraamVmxOdSZWUP969Rl4vBnZpYrhzozayhZ0iRbWmvX4sNmZjYyDnVm1pi7X83MSsOhzswaaqr71czMcuVQZ2aN1SI5snD4MzPLlUOdmTXm7lczs9JwqDOzxrq8Tp2ZmY2cQ52ZDaHrO0qYmdkIOdSZWUNN7evq7GdmliuHOjNrzGPqzMxKw6HOzBpSBMrc/er0Z2aWJ4c6M2ssgFrehTAzsywc6sysoaa2CXNDnZlZrhzqzKyx8OxXM7OycKgzs4a8TZiZWXk41JlZY26pMzMrDYc6M2usloyry8TZz8wsVw51ZjYEt9SZmZWFQ52ZNebFh83MSsOhzswa8uLDZmblMarTN5D0oKS7JS2RtDg9t7ukBZKWpT/H1r3/TEnLJT0g6ZhOl8/MhhBsmSwx7JF3Yc3MelvHQ13qXRFxSEQcmj4/A7gmIqYC16TPkXQgMB04CDgWOEfS6C6V0cwG6t9RIsthZma56laoG2gacEH6+ALghLrzcyJiQ0SsBJYDh3W/eGYGW7pfsxxmZpavboS6AK6WdJukGem5vSLiMYD0557p+QnAI3WfXZWe24qkGZIWS1q8Zs2aDhbdrMdl7np1qDMzy1s3Jkq8NSJWS9oTWCDp/iHeq0HOveRfi4iYDcwGOPTQQ/2viVmn1AL6Mv5fzP9PNDPLVcdb6iJidfrzCeAyku7UxyWNB0h/PpG+fRUwqe7jE4HVnS6jmQ1O0UQXrFOdmVmuOhrqJO0oaef+x8AfAfcA84GT07edDFyePp4PTJe0naQpwFTg1k6W0cyG4O5XM7PS6HT3617AZZL673VhRFwl6VfAXEmfAB4GTgSIiKWS5gL3ApuBUyOir8NlNLNGmglsznVmZrnqaKiLiBXA6wY5/yRwZIPPzARmdrJcZpaVW+HMzMrCO0qYWWNeg87MrDQc6sysMa9BZ2ZWGnktPmxmHfKOaf/W3gt6soSZWSm4pc6sYtTO7tJaJIeZmRWeW+rMKub6//d37btYF1vpJB0i6RZJS9IdYw6re+1MScslPSDpmJZvZmZWQW6pM7PGutut+hXgSxHxM0nHpc/fKelAYDpwELA3sFDSfl7uyMxsa26pM7PGagF9tWxH6+EvgF3Sx7uyZTeZacCciNgQESuB5SQ705iZWR2HOrOcHXHkrLyL0FjUsh+ApBlp12n/MaOJu50O/JukR4CvAmem5ycAj9S9b1V6zszM6rj71azL3nnsv3LdVV/4/fNrrzkjx9IMp7nu14iYDcxu9LqkhcArBnnpbJIFyT8dEZdI+iDwA+AoQIMXzMzM6jnUmXVZfaArvBrZZ79meFtEHNXoNUn/CZyWPp0HnJc+XgVMqnvrRLZ0zZqZWcrdr2bWWFOzX1tuPFsN/GH6+AhgWfp4PjBd0naSpgBTgVtbvZmZWdW4pc7MGuvu7NdTgG9KGgO8CMxIihBLJc0F7gU2A6d65quZ2Us51JnZELoX6iLiRuCNDV6bCczsSkHMzErKoc7MGqsF1DJuUeGpC2ZmuXKoM7PGmup+daozM8uTQ52ZDaGrY+rMzKwFDnVm1lgt2rqkiZmZdY5DnZk1FkFExjF1ZmaWK4c6M2uslu7rmomb6szM8uRQZ2aNefarmVlpONSZ2RA8UcLMrCwc6sysoagFoWwtdeGmOjOzXDnUmVljzaxT50xnZpYrhzozayyaWNLEzMxy5VBnZo1FLTmyvbmjRTEzs6E51JlZQ1FrYqycM52ZWa4c6sxsCM201JmZWZ5G5V0AMxu5I9/5zx29ftQi8+GmOjOzfDnUmZXYNded1cnL3/FkPLZlXN0Qxwu19YxydWJmlivXwmbWyH+tYTWbYuOwb3yQB5jMAV0okpmZNeJQZ2aDioiNk9iXh1k25PteiOd4lrXcyU2uT8zMcuRK2Mwauo/btxuuta6/lS7C+4mZmeXJoc6soo44clbL1xiutc6tdGZmxeGK2Kyirr3mjN8/3ufCkc+SHaq1zq10ZmbF4VBn1gM+84aFI/5so9Y6t9KZmRWLK2OzHvCpA65t6fODtda5lc7MrFgc6sxsWANb69xKZ2ZWPK6QzSyT+tY6t9KZmRWP9341s0wiYuOBeiO/5i6eZx2PssJfCs3MCsSVsplldh+3b7eWNW6lMzMrIJW9Xpa0Hngg73KMwB7A7/IuxAi43N3VznK/KiLGtelaZmZWMFXofn0gIg7NuxDNkrTY5e4el9vMzKrO3a9mZmZmFeBQZ2ZmZlYBVQh1s/MuwAi53N3lcpuZWaWVfqKEmZmZmVWjpc7MzMys5znUmZmZmVVAqUOdpGMlPSBpuaQz8i5PP0mTJP1C0n2Slko6LT2/u6QFkpalP8fWfebM9Pd4QNIx+ZUeJI2WdIekK9LnhS+3pN0kXSzp/vTP/Q9KUu5Pp/+N3CPpIkkvK0O5zcyseEob6iSNBr4LvBs4EDhJ0oH5lur3NgOfjYjXAG8GTk3LdgZwTURMBa5Jn5O+Nh04CDgWOCf9/fJyGnBf3fMylPubwFURcQDwOpLyF7rckiYAfwscGhEHA6PTchW63GZmVkylDXXAYcDyiFgRERuBOcC0nMsEQEQ8FhG3p4/XkwSMCSTluyB92wXACenjacCciNgQESuB5SS/X9dJmgj8MXBe3elCl1vSLsA7gB9AskdpRKyl4OVOjQG2lzQG2AFYTTnKbWZmBVPmUDcBeKTu+ar0XKFImgy8HlgE7BURj0ES/IA907cV6Xf5BvB5oFZ3rujl3gdYA/wo7TY+T9KOFLzcEfEo8FXgYeAx4JmIuJqCl9vMzIqpzKFOg5wr1PosknYCLgFOj4h1Q711kHNd/10kvQd4IiJuy/qRQc7l8XcwBngD8L2IeD3wHGmXZQOFKHc6Vm4aMAXYG9hR0keG+sgg5wr137yZmeWnzKFuFTCp7vlEkq6rQpC0DUmg+5+IuDQ9/bik8enr44En0vNF+V3eChwv6UGS7uwjJP03xS/3KmBVRCxKn19MEvKKXu6jgJURsSYiNgGXAm+h+OU2M7MCKnOo+xUwVdIUSduSDCCfn3OZAJAkkvFd90XE1+temg+cnD4+Gbi87vx0SdtJmgJMBW7tVnn7RcSZETExIiaT/HleGxEfofjl/i3wiKT901NHAvdS8HKTdLu+WdIO6X8zR5KMvyx6uc3MrIDG5F2AkYqIzZI+BfycZNbgDyNiac7F6vdW4KPA3ZKWpOfOAmYBcyV9guQf9BMBImKppLkkQWQzcGpE9HW91I2Vodx/A/xPGvBXAH9O8qWlsOWOiEWSLgZuT8txB8m2YDsVudxmZlZM3ibMzMzMrALK3P1qZmZmZimHOjMzM7MKcKgzMzMzqwCHOjMzM7MKcKgzMzMzq4DSLmlinSPpH4FngV2AGyJiYYP3nQD8OiLu7V7pzMzMbDBuqbOGIuKLjQJd6gTgwC4Vx8zMzIbgUGcASDpb0gOSFgL7p+fOl/SB9PEsSfdKukvSVyW9BTge+DdJSyTtK+kUSb+SdKekSyTtUHedb0m6SdKK/mumr31e0t3pZ2al5/aVdJWk2yT9UtIBXf8DMTMzKxl3vxqS3kiyLdjrSf6buB24re713YH3AgdEREjaLSLWSpoPXBERF6fvWxsR308ffxn4BPDt9DLjgbcBB5Bsd3WxpHeTtPYdHhHPp/eBZFeFT0bEMkmHA+cAR3TuT8DMzKz8HOoM4O3AZRHxPEAa1uqtA14EzpP0U+CKBtc5OA1zu5FsdfXzutd+EhE14F5Je6XnjgJ+1H/fiHhK0k4km9rPS7ZDBWC7Vn45MzOzXuBQZ/0a7heX7rN7GMmG89OBTzF4y9n5wAkRcaekjwHvrHttQ91j1f0ceN9RwNqIOKSJspuZmfU8j6kzgBuA90raXtLOwJ/Uv5i2nu0aEVcCpwOHpC+tB3aue+vOwGOStgE+nOG+VwMfrxt7t3tErANWSjoxPSdJrxvxb2ZmZtYjHOqMiLgd+DGwBLgE+OWAt+wMXCHpLuB64NPp+TnA30m6Q9K+wN8Di4AFwP0Z7nsVyfi6xZKWAJ9LX/ow8AlJdwJLgWkj/uXMzMx6hCIa9rqZmZmZWUm4pc7MzMysAhzqzMzMzCrAoc7MzMysAhzqzMzMzCrAoc7MzMysAhzqzMzMzCrAoc7MzMysAv5/nwh1qA3lo2YAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x720 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10,10))\n",
    "\n",
    "svmin = -85\n",
    "svmax = -30\n",
    "\n",
    "f = ds.frequency[1]\n",
    "\n",
    "plt.subplot(221)\n",
    "ds.Sv.sel(frequency =f ).plot(x=\"distance\", \n",
    "                          y=\"bin_depths\", yincrease=False, \n",
    "                         vmin = svmin, vmax= svmax)\n",
    "ds.depth.plot(x=\"distance\", marker='.', ls='')\n",
    "\n",
    "plt.subplot(222)\n",
    "\n",
    "ds.Sv.sel(frequency=f).where(~ds.seafloor_mask).plot(x=\"distance\", \n",
    "                         y=\"bin_depths\", yincrease=False, \n",
    "                         vmin = svmin, vmax= svmax)\n",
    "plt.title('Sv no bottom')\n",
    "\n",
    "plt.subplot(223)\n",
    "\n",
    "ds.Sv.sel(frequency=f).where(ds.swarm_mask.sel(frequency=f)\n",
    "                        ).plot(x=\"distance\", \n",
    "                               y=\"bin_depths\", yincrease=False, \n",
    "                         vmin = svmin, vmax= svmax)\n",
    "plt.title('Sv in swarms')\n",
    "\n",
    "# ds.bottom_depth.plot(marker='o', ls='')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}