{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Reading ADV Data with MHKiT" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following example demonstrates a simple workflow for analyzing Acoustic Doppler Velocimetry (ADV) data using MHKiT. MHKiT has brought the DOLfYN codebase in as an MHKiT module for working with ADV and Acoustic Doppler Current Profiler (ADCP) data.\n", " \n", "A typical ADV data workflow is broken down into\n", " 1. Review the raw data\n", " - Check timestamps\n", " - Look at velocity data quality, particularly for spiking\n", " 2. Check for spurious datapoints and remove. Replace bad datapoints using interpolation if desired\n", " 3. Rotate the data into principal flow coordinates (streamwise, cross-stream, vertical)\n", " 4. Average the data into bins, or ensembles, of a set time length (normally 5 to 10 min)\n", " 5. Calculate turbulence statistics (turbulence intensity, TKE, Reynolds stresses) of the measured flowfield\n", "\n", "Start by importing the necessary tools:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "from mhkit import dolfyn\n", "from mhkit.dolfyn.adv import api" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Read Raw Instrument Data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "DOLfYN currently only carries support for the Nortek Vector ADV. The example loaded here is a short clip of data from a test deployment to show DOLfN's capabilities.\n", "\n", "Start by reading in the raw datafile downloaded from the instrument. The `dolfyn.read` function reads the raw file and dumps the information into an xarray Dataset, which contains three groups of variables:\n", "\n", "1. Velocity, amplitude, and correlation of the Doppler velocimetry\n", "2. Measurements of the instrument's bearing and environment\n", "3. Orientation matrices DOLfYN uses for rotating through coordinate frames." ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Reading file data/dolfyn/vector_data01.VEC ...\n", " end of file at 3000000 bytes.\n" ] } ], "source": [ "ds = dolfyn.read('data/dolfyn/vector_data01.VEC')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are two ways to see what's in a Dataset. The first is to simply type the dataset's name to see the standard xarray output. To access a particular variable in a dataset, use dict-style (`ds['vel']`) or attribute-style syntax (`ds.vel`). See the [xarray docs](http://xarray.pydata.org/en/stable/getting-started-guide/quick-overview.html) for more details on how to use the xarray format." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset>\n",
                            "Dimensions:              (time: 122912, beam: 3, dir: 3, x: 3, x*: 3, earth: 3, inst: 3)\n",
                            "Coordinates:\n",
                            "  * time                 (time) datetime64[ns] 2012-06-12T12:00:02.968749 ......\n",
                            "  * beam                 (beam) int32 1 2 3\n",
                            "  * dir                  (dir) <U1 'X' 'Y' 'Z'\n",
                            "  * x                    (x) int32 1 2 3\n",
                            "  * x*                   (x*) int32 1 2 3\n",
                            "  * earth                (earth) <U1 'E' 'N' 'U'\n",
                            "  * inst                 (inst) <U1 'X' 'Y' 'Z'\n",
                            "Data variables: (12/15)\n",
                            "    beam2inst_orientmat  (x, x*) float64 2.709 -1.34 -1.364 ... -0.3438 -0.3499\n",
                            "    batt                 (time) float32 13.2 13.2 13.2 13.2 ... nan nan nan nan\n",
                            "    c_sound              (time) float32 1.493e+03 1.493e+03 ... nan nan\n",
                            "    heading              (time) float32 5.6 10.5 10.51 10.52 ... nan nan nan nan\n",
                            "    pitch                (time) float32 -31.5 -31.7 -31.69 ... nan nan nan\n",
                            "    roll                 (time) float32 0.4 4.2 4.253 4.306 ... nan nan nan nan\n",
                            "    ...                   ...\n",
                            "    vel                  (dir, time) float32 -1.002 -1.008 -0.944 ... nan nan\n",
                            "    amp                  (beam, time) uint8 104 110 111 113 108 ... 0 0 0 0 0\n",
                            "    corr                 (beam, time) uint8 97 91 97 98 90 95 95 ... 0 0 0 0 0 0\n",
                            "    orientation_down     (time) bool True True True True ... True True True True\n",
                            "    pressure             (time) float64 5.448 5.436 5.484 5.448 ... 0.0 0.0 0.0\n",
                            "    orientmat            (earth, inst, time) float32 0.0832 0.155 ... -0.7065\n",
                            "Attributes: (12/39)\n",
                            "    inst_make:                   Nortek\n",
                            "    inst_model:                  Vector\n",
                            "    inst_type:                   ADV\n",
                            "    rotate_vars:                 ['vel']\n",
                            "    n_beams:                     3\n",
                            "    profile_mode:                continuous\n",
                            "    ...                          ...\n",
                            "    recorder_size_bytes:         4074766336\n",
                            "    vel_range:                   normal\n",
                            "    firmware_version:            3.34\n",
                            "    fs:                          32.0\n",
                            "    coord_sys:                   inst\n",
                            "    has_imu:                     0
" ], "text/plain": [ "\n", "Dimensions: (time: 122912, beam: 3, dir: 3, x: 3, x*: 3, earth: 3, inst: 3)\n", "Coordinates:\n", " * time (time) datetime64[ns] 2012-06-12T12:00:02.968749 ......\n", " * beam (beam) int32 1 2 3\n", " * dir (dir) : Nortek Vector\n", " . 1.07 hours (started: Jun 12, 2012 12:00)\n", " . inst-frame\n", " . (122912 pings @ 32.0Hz)\n", " Variables:\n", " - time ('time',)\n", " - vel ('dir', 'time')\n", " - orientmat ('earth', 'inst', 'time')\n", " - heading ('time',)\n", " - pitch ('time',)\n", " - roll ('time',)\n", " - temp ('time',)\n", " - pressure ('time',)\n", " - amp ('beam', 'time')\n", " - corr ('beam', 'time')\n", " ... and others (see `.variables`)" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ds_dolfyn = ds.velds\n", "ds_dolfyn" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## QC'ing Data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "ADV velocity data tends to have spikes due to Doppler noise, and the common way to \"despike\" the data is by using the phase-space algorithm by Goring and Nikora (2002). DOLfYN integrates this function using a 2-step approach: create a logical mask where True corresponds to a spike detection, and then utilize an interpolation function to replace the spikes." ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Percent of data containing spikes: 0.73%\n" ] } ], "source": [ "# Clean the file using the Goring+Nikora method:\n", "mask = api.clean.GN2002(ds.vel, npt=5000)\n", "# Replace bad datapoints via cubic spline interpolation\n", "ds['vel'] = api.clean.clean_fill(ds['vel'], mask, npt=12, method='cubic', max_gap=None)\n", "\n", "print('Percent of data containing spikes: {0:.2f}%'.format(100*mask.mean()))\n", "\n", "# If interpolation isn't desired:\n", "ds_nan = ds.copy(deep=True)\n", "ds_nan.coords['mask'] = (('dir','time'), ~mask)\n", "ds_nan['vel'] = ds_nan['vel'].where(ds_nan['mask'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Coordinate Rotations" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that the data has been cleaned, the next step is to rotate the velocity data into true East, North, Up (ENU) coordinates.\n", "\n", "ADVs use an internal compass or magnetometer to determine magnetic ENU directions. The `set_declination` function takes the user supplied magnetic declination (which can be looked up online for specific coordinates) and adjusts the orientation matrix saved within the dataset.\n", "\n", "Instruments save vector data in the coordinate system specified in the deployment configuration file. To make the data useful, it must be rotated through coordinate systems (\"beam\"<->\"inst\"<->\"earth\"<->\"principal\"), done through the `rotate2` function. If the \"earth\" (ENU) coordinate system is specified, DOLfYN will automatically rotate the dataset through the necessary coordinate systems to get there. The `inplace` set as true will alter the input dataset \"in place\", a.k.a. it not create a new dataset." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "# First set the magnetic declination\n", "dolfyn.set_declination(ds, declin=10, inplace=True) # declination points 10 degrees East\n", "\n", "# Rotate that data from the instrument to earth frame (ENU):\n", "dolfyn.rotate2(ds, 'earth', inplace=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Once in the true ENU frame of reference, we can calculate the principal flow direction for the velocity data and rotate it into the principal frame of reference (streamwise, cross-stream, vertical). Principal flow directions are aligned with and orthogonal to the flow streamlines at the measurement location. \n", "\n", "First, the principal flow direction must be calculated through `calc_principal_heading`. As a standard for DOLfYN functions, those that begin with \"calc_*\" require the velocity data for input. This function is different from others in DOLfYN in that it requires place the output in an attribute called \"principal_heading\", as shown below.\n", "\n", "Again we use `rotate2` to change coordinate systems." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "ds.attrs['principal_heading'] = dolfyn.calc_principal_heading(ds['vel'])\n", "dolfyn.rotate2(ds, 'principal', inplace=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Averaging Data\n", "The next step in ADV analysis is to average the velocity data into time bins (ensembles) and calculate turbulence statistics. There are a couple ways to do this, and both of these methods use the same variable inputs and return identical datasets.\n", "\n", "1. Define an averaging object, create a binned dataset and calculate basic turbulence statistics. This is done by initiating an object from the `ADVBinner` class, and subsequently supplying that object with our dataset.\n", "\n", "2. Alternatively, the functional version of ADVBinner, `turbulence_statistics`.\n", "\n", "Function inputs shown here are the dataset itself; \"n_bin\", the number of elements in each bin; \"fs\", the ADV's sampling frequency in Hz; \"n_fft\", optional, the number of elements per FFT for spectral analysis; \"freq_units\", optional, either in Hz or rad/s, of the calculated spectral frequency vector.\n", "\n", "All of the variables in the returned dataset have been bin-averaged, where each average is computed using the number of elements specified in \"n_bins\". Additional variables in this dataset include the turbulent kinetic energy (TKE) vector (\"ds_binned.tke_vec\"), the Reynold's stresses (\"ds_binned.stress\"), and the power spectral densities (\"ds_binned.psd\"), calculated for each bin." ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "scrolled": true }, "outputs": [], "source": [ "binner = api.ADVBinner(n_bin=9600, fs=ds.fs, n_fft=2048)\n", "ds_binned = binner(ds, freq_units=\"Hz\")\n", "\n", "ds_binned = api.turbulence_statistics(ds, n_bin=9600, fs=ds.fs, n_fft=2048, freq_units=\"Hz\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The benefit to using `ADVBinner` is that one has access to all of the velocity and turbulence analysis functions that DOLfYN contains. If basic analysis will suffice, the `turbulence_statistics` function is the most convienent. Either option can still utilize DOLfYN's shortcuts.\n", "\n", "See the [DOLfYN API](https://dolfyn.readthedocs.io/en/latest/apidoc/dolfyn.binners.html) for the full list of functions and shortcuts. A few examples are shown below.\n", "\n", "Some things to know:\n", "- All functions operate bin-by-bin.\n", "- Some functions will fail if there are NaN's in the data stream (Notably the PSD functions)\n", "- \"Shorcuts\", as referred to in DOLfYN, are functions accessible by the xarray accessor `velds`, as shown below. The list of \"shorcuts\" available through `velds` are listed [here](https://dolfyn.readthedocs.io/en/latest/apidoc/dolfyn.shortcuts.html). Some shorcut variables require the raw dataset, some an averaged dataset.\n", "\n", "For instance, \n", "- `bin_variance` calculates the binned-variance of each variable in the raw dataset, the complementary to `bin_average`. Variables returned by this function contain a \"_var\" suffix to their name.\n", "- `cross_spectral_density` calculates the cross spectral power density between each direction of the supplied DataArray. Note that inputs specified in creating the `ADVBinner` object can be overridden or additionally specified for a particular function call.\n", "- `velds.I` is the shortcut for turbulence intensity. This particular shortcut requires a dataset created by `bin_average`, because it requires bin-averaged data to calculate.\n" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Calculate the variance of each variable in the dataset and add to the averaged dataset\n", "ds_binned = binner.bin_variance(ds, out_ds=ds_binned) \n", "\n", "# Calculate the cross power spectral density\n", "ds_binned['csd'] = binner.cross_spectral_density(ds['vel'], freq_units='Hz', n_fft_coh=2048)\n", "\n", "# Calculated the turbulence intensity (requires a binned dataset)\n", "ds_binned['TI'] = ds_binned.velds.I" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plotting can be taken care of through matplotlib. As an example, the mean spectrum in the streamwise direction is plotted here. This spectrum shows the mean energy density in the flow at a particular flow frequency." ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0.5, 1.0, 'Streamwise Direction')" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkoAAAHJCAYAAAB67xZyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB+LUlEQVR4nO3dd3wUdfoH8M9sS+8JIQmEXgy9F6UJgnoigthOEAQLJ55YsJc7/Sl6WPAKYjkB5RRQEc+zASJdpfdIFUhIAiEJ6ckmuzu/Pzazmd2dyW52N9ls8nm/XrwuOzs7803cyz55vs/3+QqiKIogIiIiIicafw+AiIiIqKlioERERESkgoESERERkQoGSkREREQqGCgRERERqWCgRERERKSCgRIRERGRCgZKRERERCoYKBERERGpYKBEFOB27tyJyZMnIzU1FUFBQUhMTMSwYcPw2GOP2Z33zjvvYPny5f4ZpJ9t3rwZgiBg8+bNjX5P6Z/BYEBCQgKuvPJKPPvsszh37pzTa5YvXw5BEHD27NlGG6dE7f1x9uxZCILQYt87RAK3MCEKXN9++y1uvPFGjB49Gvfeey+SkpKQk5ODPXv2YNWqVTh//rzt3J49eyI+Pr5Rg4Wmori4GOnp6UhLS0NkZGSj3HPz5s0YM2YMFixYgDFjxsBsNiM/Px87d+7E0qVLUVRUhA8++AB33nmn7TWXLl3C6dOn0a9fPwQFBTXKOCVq7w+j0Yj9+/ejU6dOSEhIaNQxETUFDJSIAtioUaOQlZWFY8eOQafT2T1nsVig0dQmjesTKFVXV0MQBKdrkvukQOnzzz/H1KlT7Z4rKCjAuHHjcPjwYezbtw+9evWq17XLy8sRGhrqy+G26ECaqC6ceiMKYPn5+YiPj1cMaORBUvv27XH06FFs2bLFNhXUvn17ALVTRCtWrMBjjz2GlJQUBAUF4dSpUwCAH3/8EWPHjkVkZCRCQ0Nx5ZVXYuPGjXb3OnXqFO6++2506dIFoaGhSElJwcSJE3H48GG786R7ffrpp3jyySeRlJSE8PBwTJw4ERcvXkRJSQnuu+8+xMfHIz4+HnfffTdKS0ttr7/lllvQo0cPu2tOnDgRgiDg888/tx3bt28fBEHA//73P7v7yoOA33//HbfffjuSk5NtU5Zjx47FgQMH7K6/evVqDBs2DGFhYQgPD8eECROwf/9+F/9l6hYbG4v33nsPJpMJixYtsh1XmnobPXo0evbsia1bt2L48OEIDQ3FrFmzAFgzZfPnz0eHDh1gMBiQkpKChx9+GGVlZXb3s1gs+Oc//4m+ffsiJCQE0dHRGDp0KL7++msAdb8/1Kbetm/fjrFjxyIiIgKhoaEYPnw4vv32W7tzpO9n06ZN+NOf/oT4+HjExcVhypQpyM7O9upnSNRYGCgRBbBhw4Zh586deOihh7Bz505UV1crnrd27Vp07NgR/fr1wy+//IJffvkFa9eutTvn6aefRkZGBt59913873//Q6tWrfCf//wH48ePR2RkJD766CN89tlniI2NxYQJE+yCpezsbMTFxeG1117DDz/8gMWLF0On02HIkCE4fvy403ieeeYZ5ObmYvny5XjzzTexefNm3HHHHbj55psRFRWFlStX4oknnsCKFSvwzDPP2F43btw4pKenIycnBwBgMpmwZcsWhISEYMOGDbbzfvzxR+h0OowePVr1Z3f99ddj7969WLhwITZs2IAlS5agX79+KCwstJ2zYMEC3HHHHUhLS8Nnn32GFStWoKSkBCNGjEB6enqd/21cGTRoEJKSkrB161aX5+bk5GDatGn44x//iO+++w4PPPAAysvLMWrUKHz00Ud46KGH8P333+PJJ5/E8uXLceONN0I+WTBz5kzMmzcPgwYNwurVq7Fq1SrceOONtoDMnfeH3JYtW3D11VejqKgIH374IVauXImIiAhMnDgRq1evdjr/nnvugV6vx6effoqFCxdi8+bNmDZtWv1/aET+IBJRwMrLyxOvuuoqEYAIQNTr9eLw4cPFV199VSwpKbE7t0ePHuKoUaOcrrFp0yYRgDhy5Ei742VlZWJsbKw4ceJEu+Nms1ns06ePOHjwYNVxmUwmsaqqSuzSpYv4yCOPON3L8ZoPP/ywCEB86KGH7I7fdNNNYmxsrO3xqVOnRADixx9/LIqiKG7fvl0EID7xxBNihw4dbOddc8014vDhw53uu2nTJlEUrT83AOLbb7+t+j1kZGSIOp1O/POf/2x3vKSkRGzdurV46623qr5Wfs/PP/9c9ZwhQ4aIISEhtsfLli0TAYhnzpyxHRs1apQIQNy4caPda1999VVRo9GIu3fvtjv+xRdfiADE7777ThRFUdy6dasIQHz22WfrHK/a++PMmTMiAHHZsmW2Y0OHDhVbtWpl9x4zmUxiz549xTZt2ogWi8Xu+3nggQfsrrlw4UIRgJiTk1PnmIiaAmaUiAJYXFwctm3bht27d+O1117DpEmTcOLECTz99NPo1asX8vLy3L7WzTffbPf4559/RkFBAWbMmAGTyWT7Z7FYcO2112L37t22KR6TyYQFCxYgLS0NBoMBOp0OBoMBJ0+exG+//eZ0rxtuuMHu8RVXXAEA+MMf/uB0vKCgwDb91qlTJ7Rv3x4//vgjAGDDhg3o1asXpk2bhjNnzuD06dMwGo3Yvn07xo0bp/q9xsbGolOnTnj99dfx1ltvYf/+/bBYLHbnrFu3DiaTCXfddZfd9x8cHIxRo0b5pJZHdLNENCYmBldffbXdsW+++QY9e/ZE37597cY3YcIEu2nG77//HgAwd+5cr8cLAGVlZdi5cyemTp2K8PBw23GtVovp06fj/PnzTlnEG2+80e5x7969AUBx5R9RU8NKTaJmYODAgRg4cCAAayH2k08+iUWLFmHhwoVYuHChW9dISkqye3zx4kUAcCpElisoKEBYWBgeffRRLF68GE8++SRGjRqFmJgYaDQa3HPPPaioqHB6XWxsrN1jg8FQ5/HKykrbh/LYsWPxww8/ALBOsV1zzTXo1asXEhMT8eOPP6JLly6oqKioM1ASBAEbN27ESy+9hIULF+Kxxx5DbGws7rzzTrzyyiuIiIiwff+DBg1SvIa8BsxTGRkZSE5Odnme438bwPrf59SpU9Dr9YqvkYLkS5cuQavVonXr1t4Ntsbly5chiqLimKTvJT8/3+54XFyc3WNpRZ/Se4OoqWGgRNTM6PV6/OUvf8GiRYtw5MgRt18nCILd4/j4eADAP//5TwwdOlTxNYmJiQCA//znP7jrrruwYMECu+fz8vIQHR1dj9G7NnbsWHz44YfYtWsXdu7cieeeew4AcPXVV2PDhg04d+4cwsPDVccsadeuHT788EMAwIkTJ/DZZ5/hr3/9K6qqqvDuu+/avv8vvvgC7dq18+n3AAC7du3ChQsXMHv2bJfnOv63Aaz/fUJCQrB06VLF10jjT0hIgNlsxoULFxSDm/qSgmCpTkxOKtCW7k3UHDBQIgpgOTk5ih9+0nSXPFsRFBRUr7/gr7zySkRHRyM9PR0PPvhgnecKguDU9+fbb79FVlYWOnfu7PY93TF27FgIgoDnn38eGo0GI0eOBGAt9H788cdx7tw5jBw5UjXToqRr16547rnnsGbNGuzbtw8AMGHCBOh0Opw+fdppWtJbBQUFmDNnDvR6PR555BGPrnHDDTdgwYIFiIuLQ4cOHVTPu+666/Dqq69iyZIleOmll1TPc/f9ERYWhiFDhuDLL7/EG2+8gZCQEADWlXX/+c9/0KZNG3Tt2rX+3xBRE8VAiSiATZgwAW3atMHEiRPRvXt3WCwWHDhwAG+++SbCw8Mxb94827m9evXCqlWrsHr1anTs2BHBwcF19u8JDw/HP//5T8yYMQMFBQWYOnUqWrVqhUuXLuHgwYO4dOkSlixZAsD6ob18+XJ0794dvXv3xt69e/H666+jTZs2Pv+eW7VqhZ49e2L9+vUYM2aMrZ/QuHHjUFBQgIKCArz11lt1XuPQoUN48MEHccstt6BLly4wGAz46aefcOjQITz11FMArEvmX3rpJTz77LP4/fffce211yImJgYXL17Erl27EBYWhhdffNHleE+ePIlff/0VFovF1nDyww8/RHFxMT7++GOndgfuevjhh7FmzRqMHDkSjzzyCHr37g2LxYKMjAysX78ejz32GIYMGYIRI0Zg+vTpePnll3Hx4kXccMMNCAoKwv79+xEaGoo///nPAOr3/nj11VdxzTXXYMyYMZg/fz4MBgPeeecdHDlyBCtXrlTMgBEFKgZKRAHsueeew3//+18sWrQIOTk5MBqNSEpKwrhx4/D000/biqQB4MUXX0ROTg7uvfdelJSUoF27di63ypg2bRpSU1OxcOFC3H///SgpKUGrVq3Qt29fzJw503be3//+d+j1erz66qsoLS1F//798eWXX9qmxXxNatYor0NKTU1Fly5dcPLkyTrrkwCgdevW6NSpE9555x1kZmZCEAR07NgRb775pi1wAKwtE9LS0vD3v/8dK1euhNFoROvWrTFo0CDMmTPHrbFK7Q10Oh2ioqLQtWtXzJo1C/fdd59XU3phYWHYtm0bXnvtNbz//vs4c+YMQkJCkJqainHjxtn6IAHWfkb9+/fHhx9+iOXLlyMkJARpaWl2rRfq8/4YNWoUfvrpJ/zlL3/BzJkzYbFY0KdPH3z99ddOhfpEgY6duYmIiIhUsD0AERERkQoGSkREREQqGCgRERERqWCgRERERKSCgRIRERGRCgZKRERERCrYR8kLFosF2dnZiIiIYIM1IiKiACGKIkpKSpCcnOxy30YGSl7Izs5G27Zt/T0MIiIi8kBmZqbLHQQYKHkhIiICgPUHHRkZ6efREBERkTuKi4vRtm1b2+d4XRgoeUGabouMjGSgREREFGDcKZthMTcRERGRCgZKRERERCoYKBERERGpYKBEREREpIKBEhEREZEKBkpEREREKhgoEREREalgoERERESkgoESERERkQoGSkREREQqGCgRERERqWCg5IHFixcjLS0NgwYN8vdQiIiIqAEJoiiK/h5EoCouLkZUVBSKiop8uinujlN5WLThBHqmRKFXShR6tYlCp4RwaDWuN+8jIiKiutXn81vXSGOieth37jL21PyThOi1SEuORK+UKFsA1SkhDDotk4JEREQNhYFSEzS5fwpSYkJwOKsIR7KKcDS7GOVVZuw9dxl7ZcFTsF6DtCRZ8NQmCp0Twhk8ERER+Qin3rzQUFNvjswWEWfySnE4qwiHzxfXBE9FKKsyO50brNfgiprgqReDJyIiIif1+fxmoOSFxgqUlFgsIn7PK8ORrCJrAJVVhKNZroMnadquSysGT0RE1DIxUGok/gyUlFgsIs7k1wRP52uCp+xilBpNTucG6ewzTz1TotAlMRx6Bk9ERNTMMVBqJE0tUFJisYg4m19mq3eyZp6KUaIQPBlswVNt9qlrYgSDJyIialYYKDWSQAiUlFgsIs4VlNcGT+eLcCS7CCWVKsFT6wjblJ0UPBl0DJ6IiCgwMVBqJIEaKCmxWERkFJTjkDvBk1aD7km1wVMvBk9ERBRAGCg1kuYUKCkRRRHn8svtpu2OZBWhWCV46tY6AsnRwYgJNSA61ICYUH3N13rEhFkfR4caEBWi53QeERH5DQOlRtLcAyUlomjNPB2WBU6HzysHT3WJCNIhOkyvGFRFh1gDK8fj4UE6CAK7kxMRkXcYKDWSlhgoKRFFEZkFFUjPKcKlEiMul1fjcnkVCmv+93J5NQrLq3C5rKreAZWcXisgKsQaPA1oF4O/3tgDwXqtD78TIiJqCbiFCTUqQRCQGheK1LhQl+eaLSKKKqRAqgqXy+yDqsKKasXjRpMF1WYReaVG5JUacTK3FPllVVhyZ3/2gyIiogbDQIkalVYjIDbMgNgwQ71eV1FltgVOJ3NL8PgXh7Ah/SKe/+9RLJjck1NyRETUIBgoUUAIMWgRYghBcnQI0pIjEaTT4E+f7MPKXRloHRmMeeO6+HuIRETUDHHOggLStT2T8NKkngCART+ewKpdGX4eERERNUcMlDywePFipKWlYdCgQf4eSos2fWg7PDimMwDgmbWH8WP6RT+PiIiImhuuevMCV735nyiKeOKLQ/h873kE6zX49N6h6J8a4+9hERFRE1afz29mlCigCYKABVN6YUy3BFRWWzBr+W6cyi3197CIiKiZYKBEAU+v1WDxnf3Rp200CsurMWPpLlwsrvT3sIiIqBlgoETNQqhBh6UzBqJDfBiyCiswc9luFFdW+3tYREQU4BgoUbMRFx6Ej2cNRnx4EH7LKcacFXthNJn9PSwiIgpgDJSoWWkbG4rldw9CmEGLn0/n47HPDsJi4XoFIiLyDAMlanZ6pkTh3ekDoNcK+OZQDl757jd/D4mIiAIUAyVqlkZ0ScAbt/QBAHy4/Qw+2Pq7n0dERESBiIESNVuT+qbgmeu7AwBe+e43fLU/y88jIiKiQMNAiZq1e0d0xOyrOgAAHv/iILadvOTnERERUSBhoETNmiAIePb6KzCxTzKqzSLmrNiLI1lF/h4WEREFCAZK1OxpNALeuKU3hnWMQ1mVGTOX7UZGfrm/h0VERAGAgRK1CEE6Ld67awCuSIpEXqkRM5btQn6p0d/DIiKiJo6BErUYkcF6fHT3IKREh+BMXhlmLd+N8iqTv4dFRERNGAMlalFaRQbj49mDEROqx8HzRZj7yT5Umy3+HhYRETVRDJSoxemUEI4PZw5CsF6DTccvYQEbUhIRkQoGStQi9U+Nwdu39QUArN6diSoTs0pEROSMgRK1WOPTWiM+PAjlVWbsy7js7+EQEVETxECJWiyNRsCILvEAgK0n2IiSiIicMVCiFm1kV2ugtO1knp9HQkRETREDJQ8sXrwYaWlpGDRokL+HQl66srM1UDqSXcS+SkRE5ISBkgfmzp2L9PR07N69299DIS+1igjGFUmREEVg+ylmlYiIyB4DJWrxpOm3rScYKBERkT0GStTijeySAADYdvISRFH082iIiKgpYaBELd7A9jEI1muQW2LE8Ysl/h4OERE1IQyUqMUL0mkxtGMcAGAbp9+IiEiGgRIRaqfftp5kPyUiIqrFQIkItQXdO88UoLLa7OfREBFRU8FAiQjWjXKTooJRZbJg15kCfw+HiIiaCAZKRAAEQaidfuN2JkREVIOBElGNEdzOhIiIHDBQIqpxVed4CAJw/GIJLhRV+ns4LVZRRTUy8sv9PQwiIgAMlIhsokMN6N0mGoC1+ST5x6zluzH6jU1Yd/SCv4dCRMRAiUhuVJea7Uw4/eYXoijiSFYRLCIwb9V+HDpf6O8hEVELx0CJSGZEV2tB9/aTl2CxcDuTxlZcYYLRZAEAVFZbMPujPcgurPDzqIioJWOgRCTTt200woN0uFxejSPZRf4eTotzscRaGxYRpEO3xAhcKjFi1vLdKKms9vPIiKilYqBEJKPXajC8U812Jpx+a3QXi62BUnJ0CD6cORDx4UE4dqEEy3ec9e/AiKjFYqBE5ECafmM/pcZ3sdgIAGgVGYQ2MaGYM6ojACA9p9ifwyKiFoyBEpGDUTWNJ/eeu4xSo8nPo2lZpIxSYmQwAKBTq3AAwJm8Mr+NiYhaNgZKRA5S40LRLi4UJouIX0/n+3s4LYoUKLWuCZQ6xocBsAZKLK4nIn9goESkwLadCfspNarajFIQACAlOgR6rQCjyYKcYjYBJaLGx0CJSMGILtzOxB9qa5SsGSWdVoPU2FAAwJlLnH4josbHQIlIwbBOcdBpBJzJK0NmAbfTaCy5DjVKANAh3lqn9HteqeJrDmYWYvI7O7D5eG7DD5CIWhwGSkQKIoL16J8aA4DTb43FYhGRW2LNKElTbwDQMcFap/S7Qkap1GjCgyv3YX9GIVbtymycgRJRi8JAiUiFNP3GNgGNI7+sCiaLCEEAEsJrA6UOsoJuR698+xsyC6ydu3OK2MGbiHyPgRKRipE1/ZR+PpUPk9ni59E0f1Ihd3x4EHTa2l9NHVUCpc3Hc7FyV4btcVYhi72JyPcYKBGp6JkShehQPUqMJhzk5qwNLrfEfsWbpEPN1Nv5y+UwmswAgKKKajy55hAAYHK/FABAXqkRldXmxhouEbUQDJQ8sHjxYqSlpWHQoEH+Hgo1IK1GwFWdrdNvW05w9VtDk1a8JUYE2x1PCA9CeJAOFhG2wvqv9mfhYrERHeLDsGByL4TotQCAC0XMKhGRbzFQ8sDcuXORnp6O3bt3+3so1MCkfkrbWNDd4KQgp1WkfaAkCIKtTkkq6JZWuN06sC1CDFokR1tfk13IOiUi8i0GSkR1GNHVmlE6mFmIonLuYN+QpKm31g6BEmBf0F1ZbcYvv1s7po/uZg1kk6NDAABZDJSIyMcYKBHVISkqBF1ahcMiAjtOc/qtIdmm3hxqlAD7QGnnmQJUVlvQOjIY3VtHALB28AaAbBZ0E5GPMVAicmGEtJ0J2wQ0KMcNceVsvZTyymzTbqO6JkAQBAC1GSVOvRGRrzFQInJhZNfa7UxEkRuzNpTa7UvUM0q/XyrDluPWgFWadgNkgRJ7KRGRjzFQInJhSIc4GHQaZBVW4DT3G2sQ1WYL8sukqTfnjFL7mkApr9SI3/PKoNMIuLKmISgAWzE3a5SIyNcYKBG5EGLQYnD7WABc/dZQLpUYIYqAXisgNtTg9HxksB7xsm7d/dvFIDJYb3ucIpt6Y9aPiHyJgRKRG6TtTLadZEF3Q5Dqk1pFBEOjERTPkTp0A/bTbgDQOsqaUaqstuAyVycSkQ8xUCJyg7SdyS+n823doVuiymoz1h294PMO2HXVJ0k6yAOlrq3sngvSaZEQYX0tC7qJyJcYKBG5oXvrCMSHB6Gi2oy95y77ezh+IYoi/rxyP+5fsRdLd5zx6bVt25dEONcnSaStTBIjg3BFUoTT81z5RkQNgYESkRsEQcDImum3rS10O5P/HcrBhvSLAKyZNV+SunIr9VCSXJOWiLgwA2Zd2cHWFkAuhd25iagBMFAicpM0/dYSC7rzS43469dHbY8PZBTCYvFd0XTt1Jt6RqlTQjj2Pn8N7h/VSfH5pCipRQCbThKR7zBQInLTVTUZpaPZxbhUYvTzaBrXX74+ioKyKnRvHYFQgxYlRhNO5pb67Pp1bV/iLm5jQkQNgYESkZviw4PQIzkSALDjVMuZflt39AK+OZQDrUbA61P7oE+baADA/gzf1WrV1ZXbXZx6I6KGwECJqB5a2nYmReXVeO6rIwCA+0d2RK82UeiXGg0A2OdhoCSKIv7y3yN4Y91x2/RdXfu8uYvF3ETUEBgoEdWDtJ3J1hayncn/fZuOSyVGdEoIw0NjuwAA+qfGAAD2ZxR6dM3jF0vw0S/n8K9Np/D8f4+gosqMogpr76O6apRckQKl3BIjqkwWj69DRCSn8+RFX3/9db1fc8011yAkJMST2xE1GQPaxSBEr0VeqRG/5ZQgrWYqrjnafDwXX+w9D0EAFk7tg2C9FgDQtyajdDK3FEUV1YgK0ddxFWfHckpsX3+yMwOFNQ0ig/UaRAZ79CsJABAXZoBBp0GVyYKLxZVoGxvq8bWIiCQe/Va66aab6nW+IAg4efIkOnbs6MntiJqMIJ0WQzvGYtPxS9h28lKzDpT+sfEkAGDm8PYY0C7Gdjw+PAjt4kJxLr8cBzILMaprgtolFB27YA2UOsaH4fe8Mnx7OAeAtZBbadm/uwRBQEp0CM7klSGrsIKBEhH5hMdTbxcuXIDFYnHrX2gof2FR81HbJqD5FnRXVptx6HwRAGD2VR2cnq+dfqt/ndLxC8UAgLuvbI8Xb+xhO+7NtJskmQXdRORjHgVKM2bMqNc02rRp0xAZ2Xz/8qaWRSro3nW2ABVVzXM7k0Pni2CyiEiMDLJtOCtXW9BdWO9rn7hobSvQrXUkZgxvj6ev6w5BqA2+vBEbZi0G535vROQrHk29LVu2rF7nL1myxJPbEDVJnRLCkBIdgqzCCvx8Og9jr0j095B8TtqmZUC7GMXpMCmoOZBxGRaLqLqRraPiympbn6NuidZtSO4f1Qk3D2iDuDCD1+OWapxKKhkoEZFveL3q7eqrr8aLL77odPzy5cu4+uqrvb08UZMjCALG97AGR6t2Z/p5NA1DCpTUsjzdWkcgWK9BcaUJv+e533jyRE19UlJUMKJCa4vA48ODvKpPkkQEW69ZUmny+lpERIAPAqXNmzfjX//6F2666SaUlZXZjldVVWHLli3eXp6oSbpzSCoAYONvF5tdPYwoirYeSf3bKQdKeq0GvWsaT+47V+j2taVC7m6tnTe19YWImoxScQUzSkTkGz7po/Tjjz/iwoULGDp0KM6ePeuLSxI1aZ1bRWBYxzhYRGDlrgx/D8enzuaXo6CsCgadxtaJXImtoDvT/YLu4w0cKNVOvTGjRES+4ZNAKSkpCVu2bEHv3r0xaNAgbN682ReXJWrSpg1tB8A6/dacGhzuq5l2650ShSCdVvU8W0F3PTJKUqDUvcEySjVTb0ZmlIjIN7wOlKS6gqCgIHzyySeYN28err32WrzzzjteD46oKRvfIxEJEUG4VGLE+vQL/h6Oz+zNqC3krosUKJ3ILXGreFoURRyraQ3QLbFhVsFGMKNERD7mdaDkuI3Dc889h08++QRvvvmmt5cmatL0Wg3uGNQWAPCfX8/5eTS+I2WU1OqTJK0igtEmJgSiCBzMLHJ53QvFlSiuNEGrEdCpVZhPxuooMoTF3ETkW14HSmfOnEF8fLzdsZtvvhm//vorli5d6u3liZq02wenQiMAv/5egJMXS1y/oIkrrqzG8Zrvw52+RtI57myQK+/IXdeUnjci2B6AiHzM40CpuLgYxcXFiImJQWlpqe2x9C81NRWTJ0/25ViJmpzk6BCMq+mj9MnOwC/qPpBRCFEEUmNDkRAR5PL8/jXTb+506G7oQm6gtkapuIIZJSLyDY93oIyOjq6z74koihAEAWZz8+xcTCSZNrQd1qdfxJq95/H4hG4IC/J8Y1d/kzeadEc/28q3Qtv/59U0dCE3UJtRqjJbUFlttm3kS0TkKY9/o2/atMn2tSiKuP766/Hvf/8bKSkpPhkYUaC4qnM82seF4mx+Ob4+mI07Bqf6e0gec9U/ydEVSZEI0mlQWF6NM3ll6JgQrnquNPXWNbHhAqVwgw6CAIiitU6JgRIRecvjQGnUqFF2j7VaLYYOHYqOHTt6PSiiQKLRCLhzSDu88t1vWPHLOdw+qK1Pukw3NrNFxP6avdsGuLnvmkGnQZ+20dh1pgDvbfkdf5vaW/G8arMFp3OtHby7t264fR81GgHhBh1KjCaUVFa7NX1IRFQXn/RRImrppg5ogyCdBuk5xdifWejv4XjkZG4JSo0mhBm09aojevSarhAEYPWeTPxwJEfxnLN5ZagyWxBq0KJNjPsbanuCLQKIyJcYKBH5QEyYATf0TgYQuK0CpPqkvqnR0Lq5yS0ADO0YhzmjOgEAnvryMC4UVTqdI592c3cDXU9JLQKKufKNiHzAp4FSIE43EPnK9GHWTt3fHMrB5bIqP4+m/myF3G5Ou8k9Mq4reqZEorC8GvM/PwiLxb6/2omLDV/ILWFGiYh8yeMapSlTptg9rqysxJw5cxAWZt9I7ssvv/T0FkQBpU+bKPRMicSRrGJ8vjcT943s5O8h1Yu7jSaVGHQa/P32fvjDP7Zh+6k8LN1xBveMqK1XbOjNcOVs25gwo0REPuBxoBQVFWX3eNq0aV4PhiiQCYKA6UPb4ck1h/HJzgzcc1XHBp9m8pW8UiPO5pcDqF3yX1+dEsLxwg098Mzaw1j4w3EkRAQhIlgHswU4kmXt3N04gRIzSkTkOx4HSsuWLfPlOALK4sWLsXjxYvaIIicT+yTj5W9/w7n8cmw7lYdRXRP8PSS3SNmkLq3CEVVT4+OJOwa3xabjudiQfhHzVh1wer4hV7xJpECpmIESEfmAxzVKzzzzDHbt2uXLsQSMuXPnIj09Hbt37/b3UKiJCTXocHP/NgCAFb/Ur6hbFEX87YdjuOP9X1FU0bjTRu5uhOuKIAj42829MbpbAnokR6J3myj0bRuNAe1iMG9sF8SGGXwx3Dpx6o2IfMnjjFJOTg5uuOEGaLVaTJw4EZMmTcK4ceMQFMS+JdSyTRvaDst/Poufjl1EVmEFUqJdL4cXRRH/981vWLrjDADg8z2ZdjU+Dc2b+iRHsWEGLL97sNfX8ZQto8RtTIjIBzzOKC1btgwXL17EZ599hujoaDz22GOIj4/HlClTsHz5cuTl5flynEQBo3OrcAzrGAeLCKx0c/+3RRtO2IIkAFi1OxOiKNbxCt+pMllw6Ly1hsjbjFJTEMmMEhH5kFftAQRBwIgRI7Bw4UIcO3YMu3btwtChQ/HBBx8gJSUFI0eOxBtvvIGsrCxfjZcoIEitAlbtzkSVyVLnuUs2n8Y/fjoFAHh8QjeE6LU4lVuKfTVdshtaek4xjCYLokP16Bgf5voFTRyLuYnIl3zaR+mKK67AE088gR07duD8+fOYMWMGtm3bhpUrV/ryNkRN3jVpiWgVEYS8UiPWHb2get7Hv5zF3344BgB48trumDumM67vlQQAWL3bvWyUt6T+Sf1TY5pFLzRbRsnIjBIRec+jQGnDhg0oLbXu2/TOO+/gvvvuw7Fjx+zOSUhIwOzZs/Hf//4X8+fP936kRAFEr9Xg9prNcdU6dX++JxMv/PcoAODPV3fGn0Zb+y7dPrgtAGvjylJjw2dFDtZsudKvbXSD36sxMKNERL7kUaA0f/58hIeH49dff8VHH32E0aNHY/bs2b4eG1FAu2NwW2g1AnaeKbB1ppZ8cygbT645BACYdWUHPHpNV9tzA9vFoGNCGMqrzPjmYHaDj/PQ+UIAQO9mEyhJNUoMlIjIe15NvX311Vf485//jD/+8Y8oLy/31ZiImoWkqBCMu6IVAOATWVZp428X8fCqA7CIwO2D2uL5G66wm/ISBAG3DbRmlVbtzmzQMRaWV9kaTfZpE+Xi7MBQu+qtutEK4omo+fIoUEpOTsb06dOxatUq/OEPf4DRaGTzRSIF04Zai7rX7MtCmdGEHafy8KdP9sFkETGpbzJemdxLsS5oSv820GkEHMgsxPELJU7P+4q02q19XCiiQxu+x1FjkAIlk0VEZXXdhfRERK54FCh98cUXmDx5MjZs2ICYmBgUFBTgjTfe8PXYiALelZ3i0SE+DKVGE17+Nh33fLQHVSYLxqcl4o1b+kCrssVJQkQQxl2RCABY3YBZJak+qXeb6Aa7R2MLM+gg/VjZIoCIvOVRoBQWFoYpU6agS5cuAICkpCSMHz/epwMjag40GgF3DrEWda/clYmKajNGdInHP//YD3pt3f/3u22Qdfrty/3nYTQ1TMb2YE19Up9mUp8EWH/m4UHcxoSIfMPr9gAVFRV29Unnzp3D22+/jfXr13t7aaJmYeqANgjSWf+vNrh9LN6fPhBBOq3L143smoDWkcEoLK/G+qMXfT4uURRxINM69dZc6pMk3MaEiHzF60Bp0qRJ+PjjjwEAhYWFGDJkCN58801MmjQJS5Ys8XqARIEuOtSABZN74Y9DUvHhzIEIMbgOkgBAqxFw60DrvnGf7fH99FtOUSXySo3QagT0SG5ugRJbBBCRb3gdKO3btw8jRowAYK1dSkxMxLlz5/Dxxx/jH//4h9cDJGoObh7QBgsm97JlOtx1S83qt20n85BZ4NuVpVJ9UrfECLeDt0AhNZ0sZkaJiLzkdaBUXl6OiIgIAMD69esxZcoUaDQaDB06FOfO1W/3dCKy1zY2FFd1jgdgbVDpSwdrVrz1adu8skkAM0pE5DteB0qdO3fGV199hczMTKxbt85W1J2bm4vIyEivB0jU0klF3Z/vPQ+zxXd9gaSMUp9mtOJNEhnCGiUi8g2vA6UXXngB8+fPR/v27TF48GAMGzYMgDW71K9fP68HSNTSje+RiOhQPXKKKrH15CWfXNNiEXE4S8ooRfvkmk0JM0pE5CseB0rPPPMMdu3ahalTpyIjIwN79uzBunXrbM+PHTsWixYt8skgiVqyIJ0Wk/ulAAA+81FPpd/zSlFqNCFYr0GXVuE+uWZTwkCJiHzF40ApJycHN9xwA5KSkvDCCy8gOzsb1dW1ae7Bgweje/fuPhkkUUsnTb9tSL+IvFKj19eT2gL0SomCzkU/p0AUwWJuIvIRj39DLlu2DBcvXsRnn32G6OhoPPbYY4iPj8eUKVOwfPly5OXl+XKcRC1a99aR6NM2GiaLiC/3nff6etJGuM2xPgmQ7/fGjBIRecerPyUFQcCIESOwcOFCHDt2DLt27cLQoUPxwQcfICUlBSNHjsQbb7yBrKwsX42XqMW6vSartHp3ptebvdq2LmmG9UkAG04Ske/4NOd+xRVX4IknnsCOHTtw/vx5zJgxA9u2bcPKlSt9eRuiFmlin2SEGrQ4fakMe89d9vg6RpMZ6TnFAIC+zTyjxBolIvKWrqEunJCQgNmzZ2P27NkNdQuiFiU8SIc/9ErC53vPY9XuTAxsH+vRdY7llKDaLCImVI+2sSE+HmXTIDWcLDEyo0RE3ql3RqmiokJxKu3o0aM+GRARqbt9sHX67b8HsnDyYolH15A2wu3dJhqCIPhqaE1KJDNKROQj9QqUvvjiC3Tt2hXXX389evfujZ07d9qemz59us8HR0T2+qfG4OrurVBtFvHUl4dh8aAB5QGp0WQzrU8C5DVKJq/ruYioZatXoPTyyy9j3759OHjwIJYuXYpZs2bh008/BQD+MiJqBIIg4OWbeiLMoMXec5fxn5313yboUM3WJX2b4dYlksgQa0bJbBFRavRNVmnFL2exdPsZn1yLiAJHvQKl6upqJCQkAAAGDhyIrVu34r333sNLL73UbFP4RE1NcnQInrzO2qPsb98fQ3ZhhduvLamsxulLpQCsU2/NVahBZyvovlBU6fX1iiur8cLXR/HSN+k+6WNFRIGjXoFSq1atcOjQIdvjuLg4bNiwAb/99pvdcSJqWNOGtMOAdjEoqzLjua+OuJ3RPZxVBFEEUqJDEB8e1MCj9K82MaEAgPOX3Q8k1ZzNK4P0Iz6TV+b19YgocNQrUFqxYgVatWpld8xgMGDlypXYsmWLTwdGROo0GgGvTekFg1aDn47l4n+Hctx63cFMadotugFH1zS0ibGu6Dt/udzra8mDIwZKRC1LvQKlNm3aoHXr1orPXXnllT4ZEBG5p0tiBOaO6QwAePHro7hcVuXyNbZGk22ab32SpDZQ8kVGqVz2NQMlopbEJ32UKisrcejQIeTm5sJisdg9d+ONN/riFkSk4E+jO+Hbw9k4cbEU//dtOt66tW+d59u2LmkRGSUfTr3lM6NE1FJ5HSj98MMPuOuuuxT3dhMEAWaz2dtbEJEKg06D127ujZuX/Iwv92Xhpr4pGNk1QfHc3OJKZBdVQiNYN8Nt7jj1RkS+4PUWJg8++CBuueUW5OTkwGKx2P1jkETU8PqnxmDGsPYAgGfWHkaZynL4gzVtATq3CkdYUIM15W8yfDr1Jssoncsv96h/FREFJq8DpdzcXDz66KNITEz0xXiIyAOPT+iGlOgQnL9cgbc2nAAAmMwWnMotwf8OZuONdcfx9o/W432acVsAOWnqLb+sChVVnv/RVlhehcJy61YoWo2AimozLpZ433KAiAKD139WTp06FZs3b0anTp18MR4i8kBYkA6vTO6Jmct2Y9mOM/jldD5OXSpFlcnidO5VXeL9MMLGFxWiR0SwDiWVJmQVlqNzqwiPriNNtbWODEawXoOz+eU4k1eGpCj398kzmS24XF6NhIjm3ZKBqDnyOlD617/+hVtuuQXbtm1Dr169oNfr7Z5/6KGHvL0FEblhdLdWmNwvBWv3ZyE9pxgAEGrQomtiBK5IikC3xAj0ahON/qnR/h1oI2oTE4rfcoqRebnC40BJmnZrHx+KEL0WZ/PLcTavHMPr8bfh3E/3Yd3Ri/juoRFIS470aBxE5B9eB0qffvop1q1bh5CQEGzevNmuQ7cgCAyUiBrRK5N7YkC7GCREBKF76wi0jQmFRtNyu+anRIfgt5xir+qUztS0BugQH4ZgvRY4fgln8krdfv2h84VYd/QiAODn03kMlIgCjNeB0nPPPYeXXnoJTz31FDQar0ueiMgLoQYdpg1t5+9hNBm+WPkm9U1qHxeGUIMWQG3w5I5/bDxl+5pbYhIFHq8DpaqqKtx2220MkoioyfHFyrfaqbfaQEm+Cq4up3JL8ONvF22PCytcNwUloqbF6+hmxowZWL16tS/GQkTkU75oOilllDrEh6FdbBgAILOg3K399fZlFNo9llbPEVHg8DqjZDabsXDhQqxbtw69e/d2KuZ+6623vL0FEZFHpIxSlodTb2VGE4orrX2pkqNDoKup9zKaLCgxmhAZrK/r5Th9yb6WqaiCgRJRoPE6UDp8+DD69esHADhy5Ijdc/LCbiKixta2JqOUV2rtpRRSM3XmrtwSIwDr6sHwmiadEUE6lBhNuFRidBko/X7Jmo3q0yYKB88XMVAiCkBeB0qbNm3yxTiIiHwuMkRnC2w86aWUW2xtLNlK1v8oISIIJUYT8kqM6JQQXufrf6/JKPVLjcHB80WceiMKQB7VKB06dMhp89u6HD16FCaT8rYKREQNRRAEpNRMv2V6UKckZZRaRQTbjsXXBE2XSo11vtZktiCjwDrl16+mdxUzSkSBx6NAqV+/fsjPz3f7/GHDhiEjI8OTWxEReaW+Bd05RRW4UGTNJEmBUkKkfUYJAC6V1B0oZV6uQLVZRLBeg7Qka++kwnKueiMKNB5NvYmiiOeffx6hoaFunV9VxV8OROQfthYBBa4LugvLq3D937dBq9Fgx1NjkFuiMPUW7l6gJE27tY8LQ3SoAQBQXGmC2SJC24KbgBIFGo8CpZEjR+L48eNunz9s2DCEhLi/LxIRka+0jbX+QZfpxsq3z/ecx+WaOqKM/HJcKnaeenM3oyQVcndKCEdUSG3Rd0lltS1wIqKmz6NAafPmzT4eBhFRw2hXEyidy687UDJbRHz861nb47P55bIaJYWMkosaJak1QMeEMBh0GoQZtCirMqOwnIESUSBhO20iatbaxVkDpYz8uptEbjqWi8yC2jqms3lltVNvCjVKeS4CJSmj1DHB2qRSyioVsqCbKKAwUCKiZk2aeisxmmzTako++uUsACBEX7tNidKqN3en3qRtTjrEW1sIRNVkkbjyjSiwMFAiomYtWK9F60hroHNOYY+2arMFf/vhGLadzIMgAPeN7AgAOHGxxNb3yLGPEmBtYmmxKGeoKqvNtiBLmvqLCrFWOnDlG1FgYaBERM1eqjT95rDyLb/UiFvf+wVLNp8GADwwuhNGdIkHABzMLAIAGLQaRIfWFmPHhlkzQ2aLiMsqQc/5msLx8CCd7bXRIcwoEQUiBkpE1OxJWZ2zebWBktkiYt6qA9ifUYiIYB2W3Nkfj0/ojnZx1pqiKrO1qW5CRJDddkx6rcYWLKkVdEsBWdvYUNtrpYCpiN25iQKK14HSzJkzsXXrVl+MhYioQUgF3ecKaqfe/r7xJLafykOIXosv5gzHdb2SAADx4Qbbvm6AfSG3xFUvJakovG1MbVsUFnMTBSavA6WSkhKMHz8eXbp0wYIFC5CVleWLcRER+YyUJcqoaRGw/WQe/vnTSQDAq1N6oVvr2j3gBEGwBVaAfX2SxNXKNymjlBpbe50oKaOkEChZLCLe2nAC205ecv+bIqJG4XWgtGbNGmRlZeHBBx/E559/jvbt2+O6667DF198gepq/uVERP5Xm1GyBjBLtpyCKAJ3DG6Lm/qlOJ3fPj7M9rV8xZvE1cq3TNnUm8SWUVKYeluffgH/2HgS0z/chWqz+/toElHD80mNUlxcHObNm4f9+/dj165d6Ny5M6ZPn47k5GQ88sgjOHnypC9uQ0TkkXax1sDnUokR+aVG7D5zGQBwz4iOiue3d5FRig832K6nRCmjVFvM7VwAfvpS7ZTg1hPMKhE1JT4t5s7JycH69euxfv16aLVaXH/99Th69CjS0tKwaNEiX96KiMhtUaF6W0bni73nUWW2ICU6BB1lmSO59nGyjJJSjVIdGSVRFG0b8LqbUfpdFiit3c/yBaKmxOtAqbq6GmvWrMENN9yAdu3a4fPPP8cjjzyCnJwcfPTRR1i/fj1WrFiBl156yRfjJSLyiDT99umuDADAiC7xdqvZ5NyeelOoUbpcXo1SowlA7Ya8ABBZ00eppNLk9JqTuSW2rzekX0SZ0fkcIvIPrwOlpKQk3HvvvWjXrh127dqFPXv2YM6cOYiIqC2OnDBhAqKjo729VYOYPHkyYmJiMHXqVH8PhYgaUKrDnm8juiSonivPKCUoTr1Zj+WXOk+jSfVJiZFBCK7p8g0AkcHWjFJJpX1GyWIRcfJiqe2x0WTBofNFdX8zRNRovA6U5s2bh/Pnz2Px4sXo27ev7bgoisjIsP7lFhMTgzNnznh7qwbx0EMP4eOPP/b3MIiogclXsgkCcGXnONVz48MNaBMTglCD1tasUi4uTH3Vm62HUoz96yKCrRmlsiozTLKC7azCClRUm2HQajDuikQAwJEsBkpETYXXgdJf//pXlJaWOh0vKChAhw4dvL18gxszZoxd9ouImqd2sixR7zbRiK7Ze02JIAj48oHhWPfwSFsmSE4q5i4oq4LZYRuT7EJrfVKKbNoNACJk1ymVTa1J024dE8LQt20UAOBINgMloqbC60BJbTfu0tJSBAc7z+3Xx9atWzFx4kQkJydDEAR89dVXTue888476NChA4KDgzFgwABs27bNq3sSUfPUTlZYPbJmm5K6tIoItivGloup6cxtEZ33bpOCIKl4W2LQaWwb7hZX1AZKJ2qm3Tq3CkfPFGugdJgZJaImQ+f6FGWPPvooAOtfXi+88AJCQ2t/oZjNZuzcudNuKs4TZWVl6NOnD+6++27cfPPNTs+vXr0aDz/8MN555x1ceeWVeO+993DdddchPT0dqampAIABAwbAaHROj69fvx7JyclejY+IAoc8o1RXfZI79DX7vxWWVyO/rApx4bV1TFKgFBbk/Os1IliHimozimV1SicuWjNKXRMj0CPZGiidyStDqdFk1yGciPzD4/8X7t+/H4A1o3T48GEYDLVpbIPBgD59+mD+/PleDe66667Dddddp/r8W2+9hdmzZ+Oee+4BALz99ttYt24dlixZgldffRUAsHfvXq/GIGc0Gu2CruLiYp9dm4gaVmJkEIZ2jEVltQX9UqO9vl5cmAGF5dXIKzWia2Lt9L20Yk0pyIkM0SO3xGgXKEmF3F0Tw5EQEYTWkcG4UFyJ33KKMah9rNfjJCLveBwobdq0CQBw99134+9//zsiIyN9Nih3VFVVYe/evXjqqafsjo8fPx4///xzg9zz1VdfxYsvvtgg1yaihiUIAlbdN8xn14sPD8LpS2VOK9/KjGYAQJhB6/QaqaBbahFgsYg4lWsNlLrUBFs9UyJxobgSx1QCpcyCciRGBsOg457mRI3B6/+nLVu2rNGDJADIy8uD2WxGYmKi3fHExERcuHDB7etMmDABt9xyC7777ju0adMGu3fvVj336aefRlFRke1fZmamx+MnosBW2yLAfmq/rqk3qTC8uGa/t/OXa1e8STVUCTV9my4rNKb85lA2RizchEU/nvDRd0FErniUUXr00Ufxf//3fwgLC7PVKql56623PBqYuxwbxomiqNpETsm6devcPjcoKAhBQc49VYio5YmrWfmWX+aYUVKfenPMKMlXvOm01r9bpcaUxQqb5z74qbXkYcnm03jy2u5efw9E5JpHgdL+/fttG95KtUpK6hOw1Fd8fDy0Wq1T9ig3N9cpy0RE5Gu1vZSUV70pZpRqVsJJNUrSircushonW9bJoTFlVk3bASJqXB4FSlJ9kuPXjclgMGDAgAHYsGEDJk+ebDu+YcMGTJo0yS9jIqKWw5ZRcph6K6uSAiXXNUonpRVvrcJt50jBVJFDRul/B7PtHhdXViv2eCIi3/J67WlFRQVEUbS1Bzh37hzWrl2LtLQ0jB8/3qtrl5aW4tSpU7bHZ86cwYEDBxAbG4vU1FQ8+uijmD59OgYOHIhhw4bh/fffR0ZGBubMmePVfYmIXJGaTjp257YVc7tRo3SiZurNPqMkTb3Z7/d2NNt+lW1Gfrmt7xIRNRyvA6VJkyZhypQpmDNnDgoLCzF48GAYDAbk5eXhrbfewp/+9CePr71nzx6MGTPG9liqh5oxYwaWL1+O2267Dfn5+XjppZeQk5ODnj174rvvvkO7du28/baIiOok9U5yrFGyTb0ZlAKl2oySKIo4nVsGwNps0nZOiPLUm7Q1iuRsfhkDJaJG4PWqt3379mHEiBEAgC+++AKtW7fGuXPn8PHHH+Mf//iHV9cePXo0RFF0+rd8+XLbOQ888ADOnj0Lo9GIvXv3YuTIkV7dk4jIHXFh0tRbbaBUbbagymTdx02tjxJgDYJKjSZUVFuzT8nRtbsYqNUoSZvtXpFkXWUsbe5LRA3L60CpvLzctlfa+vXrMWXKFGg0GgwdOhTnzp3zeoBERE2RlFEqNZpQWRPwlNdMuwHqnbkBa0bpUol1yi7MoEWoLPskbX0in3orqaxGQU3makTN9ivn8st89r00ZRaLiLfWH8fm47n+Hgq1UF4HSp07d8ZXX32FzMxMrFu3zlaXlJub65f+SkREjSEyWAe91rqyV5p+K60p5DZoNYoNIeXZIilQSoiwbzkitQcoqayGpWbD3cwC64q32DADeiRbf6+ebSEZpW8O5+AfP53CzGXqPe6IGpLXgdILL7yA+fPno3379hgyZAiGDbN2vl2/fj369evn9QCbosWLFyMtLQ2DBg3y91CIyE8EQbC1CJBWvpUZ1Ve8AUBETaBUUmmytRVwCpRqzrGItSvopPqktrGhiK2Z8lPqs+QPahuj+0pmQcsICKnp8jpQmjp1KjIyMrBnzx788MMPtuNjx47FokWLvL18kzR37lykp6fX2cWbiJq/2hYBNRmlOnooAfbZotySSgDOgVKwXmvLRkktAqRgoV1sKPQ1jSmrzRaffR+e2nEqD/3/bwO+O5zTYPdo6ECMyBWfbBbUunVr9OvXDxpN7eUGDx6M7t3ZOZaImi+pTinPIaOkVMgN1GaLqs0izl+2TqdJW6EonSfVKZ0rsNYjpdoFSv4PIO5ethuXy6vxwCf7GuwejJPI37xuDwAAGzduxMaNG5GbmwuLxf6vnKVLl/riFkRETU58mP02JmUuMkqhBi20GgFmi4jfL1m7cicoBUohOuSVGm0r36QapbaxITA0oYxSVSOMgXES+ZvXgdKLL76Il156CQMHDkRSUlKDbltCRNSUSFNveTWF2aV1NJsErHVNEcE6FJZX4/Qla5bIceoNcG5MKU3BxYUFQa+z/o71RaAkiiKMJguC9co1VU2BhSkl8jOvA6V3330Xy5cvx/Tp030xHiKigCEFOZecpt7UAw8pUJIKtBUDJVu/Jev1pNqn8GCdLaNkNLkOlD7bnYkz+WV4YkI3xT9iH1l9AN8fuYCtT4xBYmSwwhX8L9DjpGMXihETamiyP19yzesapaqqKgwfPtwXYyEiCiitIqwffrnFUkZJvSu3ROqTJFGuUZK2MbFmkkora2uf6lPM/cSaQ1iy+TT2ZRQqPv/VgWwYTRZ8ujPD5bX8JZDjpMyCclz79jYMWbDR30MhL3gdKN1zzz349NNPfTEWIqKA0qomGyStYHNVowQAvRy2Hak7o1QTKMmKxKUVca6KuaUeTABQWF5Vx5m1m/TWV2NUWgTyqrfDWUX+HgL5gNdTb5WVlXj//ffx448/onfv3tDr7f9aeuutt7y9BRFRk5RgC5Tc66MEAKO6JmDlrkzbY6nOSU7KOhVVWJtOyqfeNDXRidkiwmwRodVYH+85W4DzlytwU78U61iqaoMfV7FGSaVnPZkENHzGJ4DjJNt/KwpsXgdKhw4dQt++fQEAR44csXuOhd1E1JxJU28lldZtTFwVcwPAlZ3j7R4H6ZyDKnl7AHnAEx6ks5tyqzZboNVYXz/13V8AAMF6Da7tmWQLrgDX9Uzyc+tDEARbJCOKYoP8zpcXczfUPRqKJnCGSnXwOlDatGmTL8ZBRBRwIkOsU2FVJgsulRhRXlV3HyWgtju3q+sCtZvnAoBeKyBIp7Gb7qo2O69Y+2RnhjVQqrTfK64uHk+9yb4urzLXGSB6Sp5QsoiA1g/Bh8lsgU5b/0oVZpSaB580nGxpuIUJEQHWjIq8TsmdYm4AuL5X6zqflwKt0kqTXSG3IAjQyxr7VtVkiuRZpm0n82A0mVFilAdKdQdCxR5MvYmiCLMs2+NpVsr1fWq/Nlkav3fUzt/zkfbCOqz45Wy9X6tlSqlZ8EmgtG3bNkybNg3Dhg1DVlYWAGDFihXYvn27Ly7f5HALEyKSSIHSpRKjW8XcAPDazb1x68A2WHnvUMXnI2pWvZUaTbaAJ7zmmEYjQKeReilZo4gyhyDlWE6J3TGljJK8SLrUg4xSeZXZLoipb1bK3SJt+Xl+iJMwb9UBVJkteP6/R+v9WiaUmgevA6U1a9ZgwoQJCAkJwf79+2E0WosaS0pKsGDBAq8HSETUlMkLustqapTqmnoDrDVIC6f2wbBOcYrPhwdZp+fKjPKMUu2UnWOLAMdsToksEwXU9mOSM8tWxSk974pjcFafgvDjF0ow4OUf8eH2My7PlYdT/sgoedPwMhCm3vaeK8DDq/Yjq7DC30NpsrwOlF5++WW8++67+OCDD+xWvA0fPhz79jXc/j9ERE2BvJdSqRur3twhvb7EaLJdM0IWfOlrCnWqVAIleSYKUM72mCzyabP6T7053rOy2v0g5vmvjqCgrAr/9006Nh/PrfNcf2eULF6supMHShZvLtSAbl7yC746kI0rX/sJq3c33X5a/uR1oHT8+HGMHDnS6XhkZCQKCwu9vTwRUZMmr1Eqc6OY2x0RNdkjeY2SPPiq7aVksZ0nJ89EAcrZHvk+bZXVlnp/kEvZM4nRZFY501lFde25M5fVXcJg8XONkjcNEGTlZHb1XGpMZgv++MGvWPDdbx7f0xtPrjnsl/s2dV4HSklJSTh16pTT8e3bt6Njx47eXp6IqElrFVn/GiVXpHqkimqzbZ+3cNlqOWkbE6mYWymjVOYqo+TQsFI6Z83e83h3y2mXYyxxyEJVubGliqQ++9TJpwjNfsjKeNPHSZ5RcmfsW05cws+n8/H+1t89vyn5nNeB0v3334958+Zh586dEAQB2dnZ+OSTTzB//nw88MADvhgjEVGTJU29ZRVW2IqrvQ2U5Nmji8XWrt/yLJXeMaOkECjJj10ur3LK+DgGK1KG6bHPD+K174/hVG5JnWN0zChV1SP4MdUj4JGP0yyKEEURR7KKnGqkGoo3oZnd1JsbEZerbuvkH143vXjiiSdQVFSEMWPGoLKyEiNHjkRQUBDmz5+PBx980BdjJCJqsqRi7rN55bZjYQbvapSCdFoYtBpUmS3IqQmUpJVwQG0xd5VJedWbY43SsQsluP7v27DhkVHQ2FbM2Qc2Jov99Nvl8rrrlhzvaaxHjZKpjqBq64lL2JdxGQ9d3QUajWCX+TKZRWxIv4j7VuxFt8QIrHvEuezD17wp5pa3XnIno8RuAk2TT7qDvfLKK3j22WeRnp4Oi8WCtLQ0hIeH++LSRERNmlSjJGVUUmNDPWpO6Cg8WIeCsipcKFLIKDmsenOcWnOsUQKA05fKUFplsnX9dpx6M5lFu6yQqw92xyxWfTJKdWVO7lq6CwCQFBWM2walolpWl2QRRfz3QDYA4PhF9YxXeZUJxy+UoG/baK87eXtThC3YFXO7Pj8QVsm1RF79v9lisWDp0qW44YYbMHjwYNx99914+eWX8eWXXwb0RoZERO6KCw+yywTcO6KDT64rBUZKgZJBa58VclWjJKmoqp0uc8woVZstdlud3P7+r1i1S30VlOP1fV2j9PPpfAD2AZ3JItpW/NXljg92YvI7P+OLvefdHpMaX029uVPMrWmBLaB3nSnAks2nm+yqQMCLQEkURdx444245557kJWVhV69eqFHjx44d+4cZs6cicmTJ/tynERETZJWIyAu3JpVSogIwi0D2/rkurZASapRUph6kwIOKWiRd/QucRkoOWSULKJTsPPUl86roDYdy8U/N550Cs7qs+rNnRqlc/nlNefaZ7ncydYdzCwEAHy+xweBkhef3/KQzp2pNwGNm1HydXDy9cFs/PdAVr1ec+t7v+BvPxzDN4dzfDoWX/J46m358uXYunUrNm7ciDFjxtg999NPP+Gmm27Cxx9/jLvuusvrQRIRNWVtYkJwqcSI+0d2dNp7zVNS0CN9wMr7KEntAYwOq94SI4NQesm6ka5St+3yOjJKjlNvau5ebl3OL005SnydUZJaGjiuetPXY1pT9DAftOKXs4gM0WNS3xSvapTs96lzI1CSxUkNvQHwhvSLePSzAz67XnmVCQ+t3A8AGN2tFaJCXO9pKHc2r8xnY/E1jzNKK1euxDPPPOMUJAHA1VdfjaeeegqffPKJV4NrqrjXGxHJvXBDGp66rjvuGtbeZ9eUZ5AcH9dmlKwfvlKNUuso6wq8UqPZFjzdfWXtmCqqa4Mnx55EJoulXsFObonR7nFdr92XcRmvrztmyzo5fvxLpRryIm9dzTyUPOlhtoi2aceGkllQjuf/exTzVh2oGZvn15KXoLhXzC1fJef5fd1x78d7PN4MWUm1qXbA8sylu5pytY7HgdKhQ4dw7bXXqj5/3XXX4eDBg55evknjXm9EJNcvNQZzRnWyZXp8wbFppVox9+WyKpyp+Ws8MbImUKqstmVk7hiciu6tIwDYZ5SqTPafTNVm56k3d0g9nYx1vHbKOz9j8abTeHeztT+Q49SbFPBVyq4hbSgrDzZMbk69STz58L1cXmX32JuMkpxbU2+yGNAfPaO8Icj+s/jqZ9ZUePz/6oKCAiQmJqo+n5iYiMuXL3t6eSKiFs0po2Q39VazhYnJgkmLd+BodjEAoHVNoFRSabI1qowO0SOkpl2BPFByyiiZLS7rjJQW6cSGGQDUHShJdp8tqLmX/XWk+5ZX1WY4pKDBMaNUv6m3+nMMULz5yJe/1p3pxvr2XXLlxf8dxZg3NtvVkxlNZmxIv+j1tR15O/amvODP40DJbDZDp1MvcdJqtTCZGqchGBFRc+OYUUqKDrF9LQULxRXVyCio7d8kTb3llhhtAUZkiB6hNYGSfErEqT2AQjG3I6Vl/fUJlIoqqiGKzrVQ0msrq2qPS4GFxWH6Sj715irr4snqa/n9xJoGl56Sv/SaRVttX6v1kfJ1RmnZjrM4k1eGbw9l244t+PY33PvxHq+v7Ui0+7n5/PJ+5XExtyiKmDlzJoKCghSfNxqNiseJiMg1eaCUGBmkOPVW4DBNJE29SUL0WgTrtQjRW19bUFaFl/6XjvE9Ep2ClWqzRXHNldki2qbBKhUyTlKg5M60XVFFtWIAIL1WvgecdMwxoySfejOazAg1qH+MlVeZ610ULb+fKNZdK3TiYgmKKqoxqH2sy+uaLdag6+uD2Xjii0P41x/745o0+1kZ+ao3d9oJKLFYRPxw9AL6tI22HZMvMFi5K9Oj67pS38L1QOJxRmnGjBlo1aoVoqKiFP+1atWKK96IiDwkD4w6xts38JVqoS6X2QdK8eH2f7hGh1pXHkkZpUU/nsDSHWdw+/u/KjacNCpkOuTTNpXVdQRKbkwtFVVU22WlpP5TUkZJPvUmnSfPVPz420XoZBklV93Aj10owT0f1S97Il8yb3GRURq/aCtuefcXpNdMfTqzf21GQTnmrToAo8mimNWR9+Mye7idydr9WXjgk3248rWfbMeCdLWBUkMFMfLLepIMe2vDCWQXVgCwBpVr959HRn65i1c1Do8zSsuWLfPlOIiISEZeo9QxIczuOamAWr7NSGpsKHokR0Ij1H5QSUu0pUBJvspJaQsTs+iceSk1mmzXkU+NSWxTbzVB1Ib0i9BqgKu7Jzrdp7LabNdpOyxIh5JKk61GSZ5RMpqcp94+3H4GN/VNrr2eG72bNh7LdXmOnDzjJUL9Q18eQK1Pv4C05EiFc+wfP/75oTrvLdSzQaUSqVGnw0gUvvIxh8yfJ/7y9VF8cNdArN6diWfWWnt4nX3tD74YnVdaYB9QIqKmT55R6pRgn1GSulNLK7Q6JoRhy+OjEazX2r1OyiiFKOw959yZW1SsM5L3Y6pz6q1mBd69H+/BrOV7bFNn8g7eOo2Aatk9pLFK58ozVrYaJYchnZH126msx/5y7sgrNeKP/95pe6wWq5RXmTD2rS22x+42brxYUlnnfm7y4MvTZpBKs4yO05cNwdOeVXLFNQsQvpHVVDUFDJSIiJogu6k3h4yS3pZRsgZKEcF6WzYiIri20V90iDWICVUMlByLuZX7KLk99WayILuowul18iyWySLaWgPoNIKtdsZWzC0LfKoUMkqA/TSSJ+0M6vLh9jN2j9Wmqb49lIPfL9UGbGo1UI6vNlvqrpeyr/Opc6iqlK4uitZC/lvf+8Xl6/NLjVi86ZRt6xx32f+orA88CfZO5ZaqZMX8h4ESEVETJA9unDNK1l/dhWXWv8DlXbtTYmpXx8WESVNvzlUWju0B1Poo2QdKzs/HhNauerska0IpZZLKquRL02uDMZ1WQJDUYbzmuvL+SlLNk2OsIu9V5fg9qPnPr+cw5Z0dKHQofpd7d8tpLNl82q3rOX78q21m6zh2V4GD/HxPMzRKQxEh4vO9mdh1psDl6/+8cj9eX3ccd/771zrPc6zdchztO5tPoe9L63Eq17p5cbXZgqfWHHKZLfquCW5lwkCJiKgJkvcLSpa1BgBqgwVpP7ewIHlQVZt9iqrJKCltq+KUUTKLqFKYWpNPvVUoZJSkgK7KZMHF4toshBQgOW6eK/Vy0ms1tgBDqscxO+zrZraITlkdx73fHCkVXz/31RHsyyjEPzaecnpO8tr3xxSupXyuY2Dk7qI6VwkWu6k3jzNKzoORMkrukLI5py+pbylysbgSQ1/diDfXH5fdw37AC384juJKE1765jcAwGd7MrFqdyYe/HS/i/E3PQyUiIiaoN5tojBzeHv83009bcvzJQaHpothKivkHFe9yTlOo5ktFsWVa6XG2oJxpak3KWirMllwocg5o1RqtH+NdD29VmP7vqRMi+NKvGqzxSlQkme9lPo61bXhbnFltepzStSm3hw/zNXqjhyDB1cF2vJnPe3fpBa0edLQUW0M72w6hYvFRvzzp9rA037szq+5VBK4LYO8DpRmzpyJrVu3uj6RiIjcJggC/npjD0wf2s7pOb3Dfmdq9UzRIeqBkmPQUG0WFZfbl8kCHaVASZo+O36xBOcKarMQUoBU7pBRkmqW9FrBFmBYbBklh6DIbHHKrMiDOaWMUl0dsNUaPapRC1U0Dp+c7tYouWpgaTf15mlGSWEsIpQzTa6oBZ2K97CbNpQfd/8baardl7wOlEpKSjB+/Hh06dIFCxYsQFZWli/GRUREKvQOe8rZB0q1GaXImkApRGHqrbjCPoAxqWSU5MeUAin5LvFf7qv9/S9llByn66SaJ51GA01NpCQFPI4fzFUmi9MHrTyjpFSjVFeBd7XKB/9Px5S39FD7kHecenO3RsnVirNjF2r7MXkeKCmNQ/Qwo6R83DHDCdjXVCkF1O4wW0RsO5Xn0WsbkteB0po1a5CVlYUHH3wQn3/+Odq3b4/rrrsOX3zxBaqr65fmJCIi1xz3O5NPvbWVFXNLdSlKxdxKGSWlIEN+TKlGqXOrcCRFBTsdLzWasPdcAR79zH5zdKnmyaDTQFvz6S3FD04ZJZNCRkkeKClMvRWWq3/uqDVxnLVcuSmlm7XibgchrgKll7/9zfa1x8XcHr1K2Zm8MsX9/5QCJflw71+x16P77T132a2C88bmkxqluLg4zJs3D/v378euXbvQuXNnTJ8+HcnJyXjkkUdw8uRJX9ymyVi8eDHS0tIwaNAgfw+FiFogxxoleUZJp9VgYLsYGHQajOyaAEC5j5LUs0ZiUumjJP+gVMoUCIKA2welOh0vM5pw38fOH5i1GSXBlomRpt4cM0quapSUAo+8UvVamLrql5SMf3uL3eNn1h7GY58ddJp6Uq1Rcgh26pMl8mW7I1FUnx6sy4S3t2LW8t1Ox5UCavlwcxRaCwTyriY+LebOycnB+vXrsX79emi1Wlx//fU4evQo0tLSsGjRIl/eyq/mzp2L9PR07N7t/AYiImpohjqm3gBg5X1DsfuZcUiIsG5polSjJO9vBAA7Tudh+c9nnc6TPhSLKqqRVVjh9DwAhBicP0rKjCbbCjc5KVDSazW2Wh+Lwqo36d511Sgp1SPVHSjVr0bpYrH9tT7dmYE1+847bR2jWv/j2B6gPvU6bpy7L+My3t1y2i5gVJoGFCF6nGnaccq5p9EPRy7YvpaCZ7XhSsdzA7iY2+MtTCTV1dX4+uuvsWzZMqxfvx69e/fGI488gjvvvBMREREAgFWrVuFPf/oTHnnkEa8HTETU0tU19SY9HxVae45jIAU4T73JpzyevLY7yqtM+OdPp2xZphv/tR3nVPbeUmo/4LjaTVJmlBdz112jZFSoUTK6yChdKlXvlaQ0VeeJn0/b19G4O01mFt2fUHP81sqrTAjRa2ERgZO5JeiWGIEp7/wMAIgLM+CWgW0BqNUouXlTN12QtYGoqDIjWK91+TNYuSvDt4NoRF4HSklJSbBYLLjjjjuwa9cu9O3b1+mcCRMmIDo62ttbERER6l71piTFoQ8T4FzMLWfQaWC21C77N5ktdkHSiC7xKK6oxpPXdQegHCgVV1Yr1jTJM0pSrYv0Qe5YQ1RZbbZlYXQaASaLaJdFUppKc1xlJ1ffjJKadUfti79V94NzeBwepIPRpB7Iqb36t5xiXPf3bbhtYFtoNMDKXZl4ZFxX2/O/y7Z1UevM7Ukxd31GqZpRarJr2dzndaC0aNEi3HLLLQgOdi7mk8TExODMmTOqzxMRkfuc+yg5BypyGo2AqQPa4Iu9523HpIySQatxWu0W5BAoVTrUpEwb2g4TerS2PVZaVffpTuUMgjTlp9MKtroZtYxSeZXZFoTotRqYLGaXq97qqkNqsH3OXEw7Scb3aO12ZkU+1MWbrP2KVu/JtB17e+MJ29d6WZGUenuAhhX44ZA6r2uUpk+fXmeQREREvuXYHiAi2PXfvM//IQ1DO8baHksBR7De+WPAZLbY9lQzmixORdyjuyXYPZYHSkG6uj9W7DJKNZ/eZpU+ShWyjJKURZOfojSVVtf0mlKDSl/4108nseC73zxuEqlEfilXAZ7WsbFTI5K+Z7XvXRSBDIcp2y/3nVc8V8nec5c9H5yPeJ1RevTRRxWPC4KA4OBgdO7cGZMmTUJsbKzieUREVD+uapSURIXqseq+YZj7yT58K9tPK9SgQ7FDYXd2USU6xFsbVxpNFrvtL968pY/dxrSA/dRbQkQQzl9WLvoGamuU5FkxUWXVW2W12RYwOBawA8oBRF3Taw2VUSqrMuP9rb/jxj7J6JkSZTvuOO1Un01i5YXfrlbr6WqCyPxSI/7z6zmn50VRxLkC5foyd4miiPe2/o5uiRF2x6VsY10x4sjXN9k9fvSzg5jQozU+l2XI1Ny85Gecfe0P9R+wD3kdKO3fvx/79u2D2WxGt27dIIoiTp48Ca1Wi+7du+Odd97BY489hu3btyMtLc0XYyYiatHq2sLEFcdWAUoZpSuSImw9hIwmsy2jFBOqx80D2ihcs/Ya8eF1B0ryqbfa/knS/9oHOeVVZltwo1PImig1kKwrqKira7cvnMwtsQ+UHIZisohuF1aLdpkzheX4suel2banvjys+P2LAJbtOOvejVXsOJWvuB/ela/9hDuHpGLOqE6Kr1P7fnv8ZZ1X42lMXufrJk2ahHHjxiE7Oxt79+7Fvn37kJWVhWuuuQZ33HEHsrKyMHLkSK54IyLykfbxobYC7WC9BmEKDSXVONYTRQTr7R6/NqUXJvZORpC+tkZJKspWqkWyjsE+o1QXabNc66a41mNS9sRxaqyiSjb1pnOusjErBBB1bVNS3z5K9SU105QCMse7lVaqF5o7Wrz5lK02yd1xb0hX7jDuiwKirEL1jNQnOzPwaQCvanPF64zS66+/jg0bNiAyMtJ2LDIyEn/9618xfvx4zJs3Dy+88ALGjx/v7a2IiAjW6bKf5o/CV/uzEBcWpNwpWfW19sFObJjB9nV0qB63D7Y2j5SyVjvPFGD1busUidLqNsA+gIoPVw6UQvRaVFSbZXu9aWxBkNpebxWyqTfH6UZAOYDwRzG3pNos4qdjFzHnP/uwYHIvp6Bxw28qgYyCbw/l4NtDObhzSKrLtgauslTerjyzuJEJW7L5dIPcuynwOqNUVFSE3Nxcp+OXLl1CcbF135ro6GhUVbm7JJKIiFwJ0mlx26BUjEtLrNfrHIMdeaAk/zpIdt4nNSvY1AIldzJKUoaqqKYjeLBeY1uhZVFZ9SbPKDlONyqdD9RdzO2r9gB1XX/W8j2oMlkw//ODKCy3/9zT1SOglZRXmV0GG1UmC25592fV572tMc8prmwG4Y7nfDL1NmvWLKxduxbnz59HVlYW1q5di9mzZ+Omm24CAOzatQtdu3at+0JERNTgHDNKMaG1wVF8WG2QoxSYKG2FAthnlGJD9YrnSNeTPrSDdFrbXm9SbKNUo1S76s37Ym5fNZxUvb7DeOatOmD3WOl7cMWdLNgvv+dj91n11WHeftfWDumeXaWBY9NG4XWg9N5772Hs2LG4/fbb0a5dO6SmpuL222/H2LFj8e677wIAunfvjn//+99eD5aIiLwTGWIfyMSFq2WUFAIltak3WQAVHqwSKDmsWgvS1zac/L9v0rH+6AVboCE10MwpqkBBzXYhSqveHIuzf8spxtYT6rvPK7UH8OWS/k3HnGdX5HQOjUL//qP9PqiDOzivDjeaLOpbpNRwlafy9lsURfeL0B3tOtv0NrmtL68CperqakycOBGPP/448vPzbSvg8vPz8f777yMszLq8tG/fvoodu4mIqHHJgyHAPqMUKwualDJKalNv8t5JYSpZpyiHAC1Yp7XrFn3fir22jI/UF+r7IxdswY1jN3LAOdty3d+3qe5HZz3fOb3hy95KdWV1AOfNZBf9eALn8stwpqaztlIPqspqs8uu2mr/XST5dex/5w6L6N+Gks9/dcSPd/eymFuv1+PIkSMQBAHh4eHo3bu3r8ZFREQNIF4WDGkE+wAmThZEKbUNUDoG2HeDduwM3SE+DJ0SwtAvNQZHs4tl16qdepN8fTAbgPIUn9K0lTzIcSczpFTTVGlS3pOuIShtEjzq9c0AgHUPj1TM2lSZLS4DJaVsm9ybG07U+bwroijWuZqwoa349Rz+76aefru/11Nvd911Fz788ENfjIWIiBpYrKwOKTxIZ5epkWebDFrnYEVt6k2uR3LtCuh7R3TApvmj8e8Zg5z2owvSaRR3urc+53wfpQyXPEPkuA2LEqUapUqF4MUfNh3PVSzaFkXR5dSbalsAH8kuqsSL/0tv0Hs0ZV63B6iqqsK///1vbNiwAQMHDrRNt0neeustb29BREQ+YhcM6bR2mRr759wv5gaAX58ei4KyKrSNDbUdkzdfdJxWCtZroVFZBaaUuXKs7wHsM0QVbgQ8SoXeldVNo9pYrWi7ITe0ddeSzaf8OwA/8zpQOnLkCPr37w8AOHHCPr2ntDlfc7B48WIsXrwYZnPT+EuEiMhdkbJ94UwWi10AEifLNinVy9SVUWodFYzWUdZ9P1feOxT7Mi5jYu/k2uvpHQOl2oaTjoIVMkpKU2/ymh+laS1H/p56q4vJrFww3cCtn9zy6++BX5DtDa8DpU2bNrk+qZmZO3cu5s6di+LiYkRFRbl+ARFREyH/A7baZLFrVukqo+SqaFgyrFMchnWKs3+tQ/ATpNOqNspUyigpTb3Js0gV1a4DHlG0bsornwZ03PDXX8wWi2Kg5MtVeeQZ/205TEREflVtFu0KquWtApQySq6KhuuinFFSDpSU7qOUUZK2QwHcm3oDYNs0duuJS/jr10dRXOH+tiINySyKijVKFrH5zs4ECp8EStu2bcO0adMwbNgwZGVlAQBWrFiB7du3++LyRETUAKrMFrspK3mrAJ1CYFLX0ntXHAu0g/Va1UBJqS5bqT6qvJ4ZJQAoN1oDo7uW7sLyn8/inSZSf6O29crXB7Ox9cSlRh4NyXkdKK1ZswYTJkxASEgI9u/fD6PR2q+hpKQECxYs8HqARETUcORNG11ljBIjgj2+j3Mxt3qNktIKNnkQJ5EHSu7UKAFAbonRbjrr59P5br2uoantp7ayGW82Gyi8DpRefvllvPvuu/jggw+g19f24xg+fDj27dvn7eWJiMjHru/VGgDwxyGpGN2tFYZ0iMXcMZ1Uz585vD3uH9kRs65q7/E9HTNKddUoVSkUWMeEOXf8LjPWf+pt1e5M/HDkglvnNiaTpTlsH9s8eV3Mffz4cYwcOdLpeGRkJAoLC729PBER+djrU/vgxj4pGNU1AQadBqvvH6Z4XtvYEGQWVOBPozshMdLzbBKgXKOkVnvj2MEaAKJdZJSM9Vi99tWBLLfPbSwNvQ8dec7rQCkpKQmnTp1C+/bt7Y5v374dHTt29PbyRETkY2FBOlzbs7XL8zY8MgqlRhPiw4NcnuuKUnG4YxlUSnQIsgorMKlvCvZlFNo9F6Ow2W65rJhbKbhSo5bJ8jvGSk2S14HS/fffj3nz5mHp0qUQBAHZ2dn45ZdfMH/+fLzwwgu+GCMREflBsF7rdksAVxwbKkaHGuyKuR8c0xl/HtsZ+aVVyCmqdHp9dIhzRqnMWJtFqs+ebWpF5P7GybemyetA6YknnkBRURHGjBmDyspKjBw5EkFBQZg/fz4efPBBX4yRiIgCXGpsKDolhOH0pTL8Z/YQxIcH2QUseq0GQTotkqNDcKnEeRPXaIWMUkW12brFhyAo1jWpaYoZJQZJTZfXgRIAvPLKK3j22WeRnp4Oi8WCtLQ0hIeH++LSRETUDOi0Gqx7eCQ0gmDbukQesOh1guxc50BGLbNlsojQawXklVa5PRbHzXibCvaWbJp8EigBQGhoKAYOHOiryxERUTPj2JtJntiRd95Wai6p1rqg2mxBtdmCf21yvx+S0h5zaUmRSM8pdvsaviaKTWc7FbLnk0Bp48aN2LhxI3Jzc2Fx2HRw6dKlvrgFERE1M/KARR4cKU2N6TQCtBrBqdap2iQi83J5ve57ILPQ6ZjShrmO+raNxuAOsXh/6+/1up87PtnJfklNldd9lF588UWMHz8eGzduRF5eHi5fvmz3j4iISIljjZLta43zR5NGEKBTCKCMZrPiVF1dTuWWOh1zZ9WcQafB3NGd63Uv8t6kvsmuT2pAXmeU3n33XSxfvhzTp0/3xXiIiKiF0NoFSrVfaxUCH0GwBlNGh4DmaFYxkqK96/EEuBcoBek0iFIoKqeGJd+s2R+8zihVVVVh+PDhvhgLERG1IPKaankNkl4hc6QRBMXM0d3Ld+MfG0+q3kNptZwSpW1THDkWW/dtG+3WtSmweR0o3XPPPfj00099MRYiImpBtPWoUdIIgE5hSg4AvjvsvCXJn0Z3wn0jO6J76wi3xuKYqVKSXWS/KXATXTxHPub11FtlZSXef/99/Pjjj+jdu7fdfm8A8NZbb3l7CyIiaobUapQcV8dJ5+rrUYt0/8iOiA414I73f7Ud65gQht8vlSmeX2WyICpEj6KKatVrPj6+m9v3J98R4N+I1OtA6dChQ+jbty8A4MiRI95ejoiIWgj5qjf5FidKRduCoNxfSY0UbJll82U39U3BWxtOKJ4fG2bAmj8NR0ZBOW6XBVeSX56+GklRIW7fn5oPrwOlTZs2+WIcRETUwsjjobCg2oaSSgGRIAiKq+HUSMGWfEqtrkDrg7sGIjk6BMnRysGQO0FSbJgBHeLDsPccV3w3Jx7XKF1//fUoKiqyPX7llVdQWFhoe5yfn4+0tDSvBkdERM2XfNVbqKH273a1gEipEaUaKVCSr2ZTu+6ILvHomRLl9rXV7H1uHP5vUk+vr0P2/F0L5nGgtG7dOhiNtfvx/O1vf0NBQYHtsclkwvHjx70bXRO1ePFipKWlYdCgQf4eChFRwJLXKIXJAiWlztlA/abetLZAqbbbtdrrPd0k1/FVgiA0yX3kXLm+V2t/D6FJ8zhQEh3WSTo+bs7mzp2L9PR07N69299DISIKWPL4JMSgvJebnFKRt/q1awIls3zqTfn1voxtXA0xPjzIdzfzkZv6pvh7CE2a1+0BiIiIPCEPYuQ1SgBwZec4p/OV+iu5vIfd1Jvy632ZBXKVnWrsaaSEiKYXmAUajwMlQRBsEbv8GBERkTsqq2uDmGCdfaD0n9lDsP3JMXbH6rtVCWAfKKlnpJSvG6J3neVy5CroauzJlwfHuN5ypanPB/k7svB41Zsoipg5cyaCgqzRamVlJebMmYOwsDAAsKtfIiIiclRZXVs/5FiXJAgC2sSE4tuHrkJksLU/X32KuSVBOi2A6prXK3/kWhyil21PjMHR7GJsOpaL1Xsy63U/T+udGopavRe5z+NAacaMGXaPp02b5nTOXXfd5enliYiomTPKAiU1PZJrV6Mp9VeSG9Auxmlp/pJp/fHQqv147g9pMJmVcycmi/3xtrGhaBsbip9P57kcnyPX03iNm79x9TMDGj/LFWg8DpSWLVvmy3EQEVELU+nGtiFy8qmzL+YMw6zlu1FcaQIADOsYh/kTuuHmJT/bvaZfagy2PXE1AGDdUeetTgDAbFEeR4pKT6W6NLWMktZhPK9M7oln1zo2h2akVBcWcxMRkV9UupFRkjPIAqXEyGDcNqit7XGoQYsB7WLw3vQB+H7eCMXXq029VatkmmYMb4/J/VLwzzv6KT6vVJdbj56YAIB3pw1o0IJrxwzXlH5tGuxegP97HjUEBkpEROQX49Os/XuSo4LdOl9ezC0I9sXWUnuBCT1a44qkSOXXq0QxZotyoBSs12LRbX0xsU+yW+MDnDM4ANBL1sxSPs0lCMC1PVvbbd/ia47fm1Ig09Sn3vwdfHm9hQkREZEnruwch2/+fBVS40LdOl8e6AiCgGBZ7yV3VqipFYOrBUqecLnqTfa1lCFzp47IVzSCgB7JkTiaXdwg1xfQ/CbymFEiIiK/EAQBPVOibKvaXJEHFBoBiJC9LtSNhpVqU28mlRolV+SNlge0i7GOSyHoiQzRYUy3BADA9KHtal9f878NWdfkGLgJAtC5VbjdseYW2PgaAyUiIgoIWq08UBIQHVIbKIUYXE+QqPVRUlsN5453p/VHp4QwPDC6k3WMCkFPq4hgvHPnAKy8dyj+fLWsr1HNbR3bE/iS42yjgIadamuO/RQ59UZERAFBHoQIAhAlD5TcmHpTm+LyZurt2p5JuLZnUu0YFe4RrNcixKDFsE723cbFmkjJ3ECRy7grEiE4tGvUCAIiQ+w/+pt+jZJ/gy9mlIiIKCDIgxABAqJD6zv1pvyRd7m8yvvB1VCaRlP7nJcCFE9m/r58YLjLc6JD9U73FgTgkXFd7cfRxCff/L2XLAMlIiIKCFqHGiV5Rkn+tRq1LVDySn0ZKDkfU8uHSB//9c1o9UyJRP/UGLSOVF4tOKZbAtrHhWL++G7OYxEExDXBjXnr4u/u4px6IyKigCD/vNQIAuJlH/hjurdy+Xq9SnuA7q0j6jWOuDAD8suqMPaKRKfnlKbe1DJKUm1SfafepOJ3eSZoYLsY7KnpSn7vyI4Y3inerWtN6ZfS5KfeHPcBbGwMlIiIKCDIa1U0goCwIB1W3zcUOq3GraaNWoeM0uD2sUhLjsTdV7av1zi+nzcCO88U4Nqeresco+2YSk6pdurNdaQSrNfYbSIM2BehL5jSC+MXbXW6n6v6nqhQ91Yc+pO/u51z6o2IiAKC3edlzddDOsbZlua74hiQJEUH46839kC7uLB6jaNVZDAm9kl2e5NeV5/z5VV1dyifPrQd/vfgVbbHUoBVZa4NnBxbJ7hLFJt+ewB/11AxUCIiooCgETwLBiQRwf6ZRHE11AqFrVzkWa7/u6knuiTWTg9KgUO1LFASBOUskssten0w79ZO1jC0rvvVZ++8P/SqXUnY0NuuuMJAiYiIAoJjjVJ9RYca8NGswT4ckXscp7+kgE2qjeqXGm33/Kf3DkGP5Ci4It+jTn4H+c/J1Y/JF7kad1ozAMB/H7zS7WuO75GI3xdcj0N/He925/aGwkCJiIgCgrz2xtO6lVFdE3w1HI99+afhuHVgG3xw10AAQJjBua9RXd+dlATSCso/j/r8aETRs6ySoBK01nXv+v4X02gEt7u2NyQGSkREFBDqkylxR2Ot9nIca5fECCyc2gdtY0MVn7eIolvf3/JZg9AqIggf3DXQ4Xz7flN1ESFiRJf6B48vTepp+9rV/nYAEKTT1Ktx5OAOsfUeU0NhoERERAFB8DBr4m+ughVHZouI2DCD6vNSfDe8Uzx2PTsO16TZtymo19SbCMSGGXDkxQn1GmOwrjZ8cCdQEgT3M0q/Pj0WSVHu1zM1NLYHICKigKBRmWpq6lwN1fF7sYgixnRrhdlXdcAVSZHOL1DIhAl2wZH7xdzSQsDwoPqFA/Ixq20N46nWUcqNNP2FgRIREQUEtboYTzXWonNXI50/vhu2nLhke2yxWIOd529Ic/seaisCXf2YeqYoBGJuiJf1rXLnP4VGEAIqCyjHqTciIgoIdgGAD66nsqOJz7kKEHq1icJvL11re2xxUTyl1FfILqPkxk/n+3kj8NeJabh9UKrLc5UM61i7wa/JjYaZS2cOqvcUZFPBQImIiAKCr2qUnri2G1pHBmP+BOe90BqCO0XMIbJNfV0FSor3gNrPRvneVyRFYuaVHdyqL1K8n+xl8r3qlIIhg06DoR3jfBPd+gEDJQ8sXrwYaWlpGDRokL+HQkTUYtgvQ/f8U/eB0Z3xy9NXo01M4/Tnqe9IXSVolOIooR7TbfWlVOgtv8Udg61ZqbHdW6F9vPPPNKim8DtQp95Yo+SBuXPnYu7cuSguLkZUlOumYERE5D1fftB6E2jV1wSFPeHq4nrqzZna1Ju33+bfbu7lVOjdMSHM7ufXLjYUR16cgDCDFufyy/Ha98cwuX8K7l+xFwAQXNOQMkDjJAZKREQUGHy8uKrBTRuaipnDO6Bzq/B6vc7sRs2PI7tmnBr5ce/cOrCt0zHHVW6CINiCqfbxYXh3+gC7JpYR9VxR19Rw6o2IiAJCILUEAKwdt+sbJAHqjTDvHGKd4nrsmq5Oz6lnlNz/mfVItl8Blxobqvh6jWBfiaR0C/nruidF1HssTQkDJSIiCgiB9kHrafsBtam3l2/qiYN/GY/hneOdnlMLIuXZH1ebAr83fQBuHVi7Aa3S6jrA2mDSnVYN3z50FWYOb4+Xb+oFoO7s1tPXda9zbP7EQImIiAJCYIVJQP/UGI9epzb1JggCokKU9z5T+9lc1SUevdtEYWC7GOx46uo679smJhQLp/axPVbLbFkDJeW+TXI9kqPw1xt72LqMO8ZTI7pYA77U2FBEh/p/Tzc1gT1xSERELUag1Chtnj8aR7OLMaFHouuTFXjUHkBl1Zteq8HXD17l0TjUhuGYQXI30+fYOuD1qX3w3wNZuLFvMrbKGm42NQyUiIgoIGgCJFJqHx+G9vFhHr8+Qdb12l2N2cxx2tB29vd289aO58WE6XH/qE4+GlXD4dQbEREFhECrUaqv96cPwENXd8aYbq3q/VrBh5/mD43tAgB4aVIPp+emDmiDm/un2B3ztMg+UIrzmVEiIqKAEBgfq54b36M1xveoX88liS9/No9e0xVzRnVEqME5RLgiKdIpYHU30ec4pSgPlDyYbWw0zCgREVFACJQMhD/4OtumFCQ5ahcXihC9Fl0TI9y6pk5jH3IEyEwqM0pERBQYAuWD1R8Ela8b8j4bHx0FsygiSKdVPV/OoLMPlHy1d19DY6BEREQBgRkldf742ei0mhYRRHDqjYiIAgPjJFWNFSe1xFiVgRIREQUExy02iBpDS8iaERFRM9AjOQofzxqM5OgQfw+lyVFrOOnz+zTcpZssBkpERBQwRnZN8PcQmqRAWWofiDj1RkREFOAaK9PT3Jt+KmGgREREFOACZal9IGKgREREFOACpcfUzf3b+HsI9cYaJSIiogAnCAKuSUvE5bIqdIwPb8D7ePf6YH3g5WcYKBERETUDH9w10N9DcCkQpwUDL7QjIiIivwjAOMdrDJSIiIjIr7q3brrNRDn1RkRERO7xcu5MUMlJ9WkbjQ9nDETb2FCvrt8QGCgRERGR3429ItHfQ1DEqTciIiJyC2uUiIiIiFR4u2qNq96IiIiImhEGSkREROQWtWLs5oyBEhEREZEKBkpERETklvZxTW/5fkNjewAiIiKq02f3D8PxiyUY1inOq+sE4sQdAyUiIiKq0+AOsRjcIdbfw/ALTr0RERERqWCgRERERI1CCMBGSgyUiIiIiFQwUCIiIiJSwUCJiIiISAUDJQ8sXrwYaWlpGDRokL+HQkRERA2IgZIH5s6di/T0dOzevdvfQyEiIqIGxECJiIiIGkUALnpjoERERESkhoESERERkQoGSkRERNQohADc7Y2BEhERETUK1igRERERNSMMlIiIiIhUMFAiIiIiUsFAiYiIiBpFAJYoMVAiIiIiUsNAiYiIiEgFAyUiIiJqFGwPQERERNSMMFAiIiIiUsFAiYiIiBqFEIBzbwyUiIiIqFGM6poAANBqAidg0vl7AERERNQyXNk5Hl/MGYb28WH+HorbGCgRERFRoxnYPtbfQ6gXTr0RERERqWCgRERERKSCgRIRERGRCgZKRERERCoYKBERERGpYKBEREREpIKBEhEREZEKBkpEREREKhgoEREREalgoERERESkgoESERERkQoGSkREREQqGCgRERERqdD5ewCBTBRFAEBxcbGfR0JERETukj63pc/xujBQ8kJJSQkAoG3btn4eCREREdVXSUkJoqKi6jxHEN0Jp0iRxWJBdnY2IiIiIAiCv4fjZNCgQdi9e3dA3cfTa9X3de6e7855rs5Rer64uBht27ZFZmYmIiMj3Rt0E9BY7ylf36u5va/UnuP7qvHu5c11/PW+8vT55vi+EkURJSUlSE5OhkZTdxUSM0pe0Gg0aNOmjb+HoUqr1TbKm9qX9/H0WvV9nbvnu3Oeq3Pqej4yMjKgfvE01nvK1/dqbu8rV6/n+6rh7+XNdfz1vvL2+eb2vnKVSZKwmLsZmzt3bsDdx9Nr1fd17p7vznmuzmms/w6NoTG/F76v1M9pTu8pIDDfV95cx1/vq5b0uwrw3ffDqTciPyguLkZUVBSKiooC6i80atr4vqKG0NLfV8woEflBUFAQ/vKXvyAoKMjfQ6FmhO8raggt/X3FjBIRERGRCmaUiIiIiFQwUCIiIiJSwUCJiIiISAUDJSIiIiIVDJSIiIiIVDBQIgoAkydPRkxMDKZOnervoVAA++abb9CtWzd06dIF//73v/09HGoGWsLvJrYHIAoAmzZtQmlpKT766CN88cUX/h4OBSCTyYS0tDRs2rQJkZGR6N+/P3bu3InY2Fh/D40CWEv43cSMElEAGDNmDCIiIvw9DApgu3btQo8ePZCSkoKIiAhcf/31WLdunb+HRQGuJfxuYqBE5KWtW7di4sSJSE5OhiAI+Oqrr5zOeeedd9ChQwcEBwdjwIAB2LZtW+MPlAKat++z7OxspKSk2B63adMGWVlZjTF0aqL4u8s9DJSIvFRWVoY+ffrgX//6l+Lzq1evxsMPP4xnn30W+/fvx4gRI3DdddchIyPDds6AAQPQs2dPp3/Z2dmN9W1QE+ft+0ypykIQhAYdMzVtvvjd1SKIROQzAMS1a9faHRs8eLA4Z84cu2Pdu3cXn3rqqXpde9OmTeLNN9/s7RCpGfDkfbZjxw7xpptusj330EMPiZ988kmDj5UCgze/u5r77yZmlIgaUFVVFfbu3Yvx48fbHR8/fjx+/vlnP42Kmht33meDBw/GkSNHkJWVhZKSEnz33XeYMGGCP4ZLAYC/u2rp/D0AouYsLy8PZrMZiYmJdscTExNx4cIFt68zYcIE7Nu3D2VlZWjTpg3Wrl2LQYMG+Xq4FKDceZ/pdDq8+eabGDNmDCwWC5544gnExcX5Y7gUANz93dUSfjcxUCJqBI61IKIo1qs+hKuTyB2u3mc33ngjbrzxxsYeFgUwV++plvC7iVNvRA0oPj4eWq3WKXuUm5vr9Jcakaf4PiNf43uqFgMlogZkMBgwYMAAbNiwwe74hg0bMHz4cD+Nipobvs/I1/ieqsWpNyIvlZaW4tSpU7bHZ86cwYEDBxAbG4vU1FQ8+uijmD59OgYOHIhhw4bh/fffR0ZGBubMmePHUVOg4fuMfI3vKTf5edUdUcDbtGmTCMDp34wZM2znLF68WGzXrp1oMBjE/v37i1u2bPHfgCkg8X1Gvsb3lHu41xsRERGRCtYoEREREalgoERERESkgoESERERkQoGSkREREQqGCgRERERqWCgRERERKSCgRIRERGRCgZKRERERCoYKBERERGpYKBERORHM2fOhCAIEAQBX331lU+vvXnzZtu1b7rpJp9em6ilYKBERD4l/+CX/5Nvvkn2rr32WuTk5OC6666zHVMLnGbOnOl20DN8+HDk5OTg1ltv9dFIiVoenb8HQETNz7XXXotly5bZHUtISHA6r6qqCgaDobGG1WQFBQWhdevWPr+uwWBA69atERISAqPR6PPrE7UEzCgRkc9JH/zyf1qtFqNHj8aDDz6IRx99FPHx8bjmmmsAAOnp6bj++usRHh6OxMRETJ8+HXl5ebbrlZWV4a677kJ4eDiSkpLw5ptvYvTo0Xj44Ydt5yhlYKKjo7F8+XLb46ysLNx2222IiYlBXFwcJk2ahLNnz9qel7I1b7zxBpKSkhAXF4e5c+eiurrado7RaMQTTzyBtm3bIigoCF26dMGHH34IURTRuXNnvPHGG3ZjOHLkCDQaDU6fPu39D9bB2bNnFbN3o0eP9vm9iFoqBkpE1Kg++ugj6HQ67NixA++99x5ycnIwatQo9O3bF3v27MEPP/yAixcv2k0XPf7449i0aRPWrl2L9evXY/Pmzdi7d2+97lteXo4xY8YgPDwcW7duxfbt2xEeHo5rr70WVVVVtvM2bdqE06dPY9OmTfjoo4+wfPlyu2DrrrvuwqpVq/CPf/wDv/32G959912Eh4dDEATMmjXLKZO2dOlSjBgxAp06dfLsB1aHtm3bIicnx/Zv//79iIuLw8iRI31+L6IWSyQi8qEZM2aIWq1WDAsLs/2bOnWqKIqiOGrUKLFv37525z///PPi+PHj7Y5lZmaKAMTjx4+LJSUlosFgEFetWmV7Pj8/XwwJCRHnzZtnOwZAXLt2rd11oqKixGXLlomiKIoffvih2K1bN9FisdieNxqNYkhIiLhu3Trb2Nu1ayeaTCbbObfccot42223iaIoisePHxcBiBs2bFD83rOzs0WtVivu3LlTFEVRrKqqEhMSEsTly5fX+fOaNGmS03EAYnBwsN3PMSwsTNTpdIrnV1RUiEOGDBFvuOEG0Ww2u3UPInKNNUpE5HNjxozBkiVLbI/DwsJsXw8cONDu3L1792LTpk0IDw93us7p06dRUVGBqqoqDBs2zHY8NjYW3bp1q9eY9u7di1OnTiEiIsLueGVlpd20WI8ePaDVam2Pk5KScPjwYQDAgQMHoNVqMWrUKMV7JCUl4Q9/+AOWLl2KwYMH45tvvkFlZSVuueWWeo1VsmjRIowbN87u2JNPPgmz2ex07uzZs1FSUoINGzZAo+FkAZGvMFAiIp8LCwtD586dVZ+Ts1gsmDhxIv72t785nZuUlISTJ0+6dU9BECCKot0xeW2RxWLBgAED8Mknnzi9Vl5ortfrna5rsVgAACEhIS7Hcc8992D69OlYtGgRli1bhttuuw2hoaFufQ+OWrdu7fRzjIiIQGFhod2xl19+GT/88AN27drlFAgSkXcYKBGRX/Xv3x9r1qxB+/btodM5/0rq3Lkz9Ho9fv31V6SmpgIALl++jBMnTthldhISEpCTk2N7fPLkSZSXl9vdZ/Xq1WjVqhUiIyM9GmuvXr1gsViwZcsWp0yP5Prrr0dYWBiWLFmC77//Hlu3bvXoXu5as2YNXnrpJXz//fcNUgdF1NIxP0tEfjV37lwUFBTgjjvuwK5du/D7779j/fr1mDVrFsxmM8LDwzF79mw8/vjj2LhxI44cOYKZM2c6TS9dffXV+Ne//oV9+/Zhz549mDNnjl126M4770R8fDwmTZqEbdu24cyZM9iyZQvmzZuH8+fPuzXW9u3bY8aMGZg1axa++uornDlzBps3b8Znn31mO0er1WLmzJl4+umn0blzZ7spQ187cuQI7rrrLjz55JPo0aMHLly4gAsXLqCgoKDB7knU0jBQIiK/Sk5Oxo4dO2A2mzFhwgT07NkT8+bNQ1RUlC0Yev311zFy5EjceOONGDduHK666ioMGDDA7jpvvvkm2rZti5EjR+KPf/wj5s+fbzflFRoaiq1btyI1NRVTpkzBFVdcgVmzZqGioqJeGaYlS5Zg6tSpeOCBB9C9e3fce++9KCsrsztn9uzZqKqqwqxZs7z4ybi2Z88elJeX4+WXX0ZSUpLt35QpUxr0vkQtiSA6TuoTEQWA0aNHo2/fvnj77bf9PRQnO3bswOjRo3H+/HkkJibWee7MmTNRWFjo8+1LGvseRM0VM0pERD5iNBpx6tQpPP/887j11ltdBkmSb775BuHh4fjmm298Op5t27YhPDxcsYCdiNzDYm4iIh9ZuXIlZs+ejb59+2LFihVuvWbhwoV47rnnAFhX+fnSwIEDceDAAQBQbL9ARK5x6o2IiIhIBafeiIiIiFQwUCIiIiJSwUCJiIiISAUDJSIiIiIVDJSIiIiIVDBQIiIiIlLBQImIiIhIBQMlIiIiIhX/D4kxtv/JZ6FUAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "%matplotlib inline\n", "\n", "plt.figure()\n", "plt.loglog(ds_binned['freq'], ds_binned['psd'].sel(S='Sxx').mean(dim='time'))\n", "plt.xlabel('Frequency [Hz]')\n", "plt.ylabel('Energy Density $\\mathrm{[m^2/s^s/Hz]}$')\n", "plt.title('Streamwise Direction')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Saving and Loading DOLfYN datasets\n", "Datasets can be saved and reloaded using the `save` and `load` functions. Xarray is saved natively in netCDF format, hence the \".nc\" extension.\n", "\n", "Note: DOLfYN datasets cannot be saved using xarray's native `ds.to_netcdf`; however, DOLfYN datasets can be opened using `xarray.open_dataset`." ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "# Uncomment these lines to save and load to your current working directory\n", "#dolfyn.save(ds, 'your_data.nc')\n", "#ds_saved = dolfyn.load('your_data.nc')" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3.9.12 ('base')", "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.12" }, "vscode": { "interpreter": { "hash": "357206ab7e4935423e95e994af80e27e7e6c0672abcebb9d86ab743298213348" } } }, "nbformat": 4, "nbformat_minor": 4 }