{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Finite-Time Lyapunov Exponents\n",
    "\n",
    "Code taken from <a href=\"https://github.com/LauraGomezNavarro/OceanParcels_Lyapunov\" target=\"_blank\">here</a> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wed Jan 31 13:59:40 2024\n"
     ]
    }
   ],
   "source": [
    "import time\n",
    "print(time.ctime(time.time()))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 0. Import packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import xarray as xr\n",
    "import matplotlib.pyplot as plt\n",
    "import cartopy.crs as ccrs\n",
    "import cartopy.feature as features\n",
    "from matplotlib import colors\n",
    "import sys\n",
    "from tqdm import tqdm"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Loading data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "home_folder = '/nethome/manra003/Lagrangian_diags/'\n",
    "output_folder = home_folder + 'outputs/ftle/'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Td: Time in days selected to get the location of particles after release**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "Td=15"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.path.insert(0, home_folder + \"Diagnostics/Functions/\")\n",
    "import FTLE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "ds = xr.open_dataset(home_folder + 'Simulations/toy_data_01.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get release particles from the meshgrid for release particles "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# hard coded here, as its not possible to tell from the parcels output files itself.\n",
    "coords0, coords1 = 12, 12\n",
    "\n",
    "# initial position\n",
    "x0 = np.reshape(ds['lon'][:,0].data, (coords0, coords1))\n",
    "y0 = np.reshape(ds['lat'][:,0].data, (coords0, coords1)) \n",
    "\n",
    "# final position\n",
    "x1 = np.reshape(ds['lon'][:,-1].data, (coords0, coords1))\n",
    "y1 = np.reshape(ds['lat'][:,-1].data, (coords0, coords1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "H = x0.shape[0]\n",
    "L = x1.shape[1]\n",
    "FTLE_f = np.ones_like(np.asarray(x0))\n",
    "FTLE_f[:,:] = np.NaN"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Estimate FTLE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.13078183 0.27300894\n"
     ]
    }
   ],
   "source": [
    "ftle_array = FTLE.compute_ftle(x0, y0, x1, y1, Td)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Plotting\n",
    "Plot FTLE values at the release locations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm4AAAFICAYAAADptXKlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABd00lEQVR4nO3deVxU9f7H8deZgWFRwAVBNMKlcm1xKffKXNJrlm1qpmZhZvnLpbqVmUt21TbLvKVXLTMtkSy1jVLLMk3LMi3LXVPUQAUVkG1g5vz+GBwjsAAHh+X9vI/zaPjyPd/zmbkP8MN3NUzTNBERERGRMs/i7QBEREREpGiUuImIiIiUE0rcRERERMoJJW4iIiIi5YQSNxEREZFyQombiIiISDmhxE1ERESknFDiJiIiIlJO+Hg7ABEREamYsrKysNvtJbrXZrPh7+/v4YjKPyVuIiIi4nFZWVnUDAgkg5Id0FS7dm1+//13JW9/ocRNREREPM5ut5OBySCqYMMo3r2YLEpMxG63K3H7CyVuIiIiUmr8sWAzipe4WXSM+jlpcYKIiIhIOaEeNxERESk1ForfS6RepXNT4iYiIiKlxjDAUryRUteMOI2WFqrIidv5LOkVERER7/D2thrqcfOsIiVuWVlZ1K9fn8TExNKOR0RERDzI29tqWAwDS3EXJ4B63M6hSImb3W4nMTGRQ4cOERwcXNoxlTlOp5OaNWvy7LNPMHDgnd4OR0RKQW5uLl9+uY6YmGV8/fUGAgICuP322xkyZAitWrXCKOY/PFK5nDx5kpdffpmZM2cyYcKjREffXWrPOn06nR9+2MrGjT+wadMWUlLScDgcOBwOjh49Tk5ODldf3ZrHH3+CNm3acPHFF3t1Ww31uHlWsea4BQcHV8rELTk5GafTyUUX1SE4OMjb4YiIB2VlZfP66/OZN28RR48ep3Xr1syZM4f+/fsTFKSfdyma4OBgXn31VWw2G5Mnv8TatRvp2PEaoqIiiYqKpFGjhoSEnN+/n/HxR3j99TdZtGgpmZlZ1KlTh+uvv546depgsViwWq1ERETQrl079x8bqampHnqHUlZocUIRmHn7yVitVi9HIiKe9Pnnaxg7dgqHD/9BdHQ0w4cP56qrrvJ2WFKOPf/881x++eW8884i/ve/RSQnJwNgsVho06YVjz8+guuv71CsNg8ciGfatJl88MEnhISE8PjjTzBgwAAuvfTSctETbCnB4gT1uJ2bErciCAoKwmq1cvTocW+HIiIesH//QcaOncKqVV/RrVtXPvtsJY0bN/Z2WFIBWCwWBg8ezODBgwFIS0vj999/5/vvv+ftt9/m1luHMGzYICZO/DeBgQF/29bJk6eYPn02c+cuolatWrz88stER0dTpUqVC/FWPEZDpZ6lz6YI/Pz8uOqqK/nhh63eDkVEzoPT6eSFF/5Lu3b/YufOfSxbtoyVK1cpaZNSExQUxBVXXMH999/PN998w6uvvsrChUvp3PlWfvrpl0LvMU2T2NgPadWqG2+//R4TJ05kz549jBw5stwlbQCGYZToksKpx62I2rVrz8qVcd4OQ0RKKCcnh//7v7EsXfoRTz31FE899RSBgYHeDksqEYvFwsiRI+nWrRuDBg2ie/e+/OtfXbj00obY7XbsdjvZ2XY2bdrCjh27GTBgANOnT6d27dreDv28qMfNs5S4FZFpmvj6+no7DBEpgZycHAYMeJC1azewZMkS+vbt6+2QpBJr0qQJGzdu5PXXX+e992L59dfPsNls7uuqq1oxc+brdO3a1duheoTmuHmWErciSkhIoHbtWt4OQ0RK4OOPV/HFF2v5/PPPufHGG70djgi+vr6MHj2a0aNHezsUKWeUuBWRxWLh1CktqxYpjxYufI+OHTsoaRPxAoPi96Bphtu5qTeyiO68805++ukX9uzZ7+1QRKQY1qxZx9q1G7j//mHeDkWkUjpzckJxLymcErciuvnmmwkJCWHJkhXeDkVEiigx8RjR0aP51796cvfdpbeTvYicm6WElxROn00R+fv706dPH1avXuvtUESkiJ5+eho2mz+LFr2jDbRFvOTM4oTiXlI4JW7FcMkll5CQcNTbYYhIMZimicPhAGDHjh2MGDGC33//3ctRiVQe6nHzLH02xVCnTh2SkpKx2+3eDkVEimDKlKcwTSe33XYbx48fZ/DgwcyaNYv169d7OzQRkRJR4lYMderUAeDo0SQvRyIiRREeXot3353N9u2/ERYWxo8//ghAvXr1vBuYSCViwSjRJYVT4lYMZxK3w4f/8HIkIlJU11zTgvXrPyYsLNRddvHFF3sxIpHKRXPcPEuJWzE0bNiQ2rVr89BDT7Bz5x5vhyMiRRQREc62bWuZOPExAGrV0mbaIhfKhZzjNmvWLOrXr4+/vz+tWrVi3bp156y7bNkyunXrRq1atQgODqZdu3asXLmyQL1Tp04xYsQIIiIi8Pf3p0mTJsTFee8ITCVuxVClShU2btxI1apB9Ox5F8nJJ7wdkogUkc1mw9/fj4CAAJ1RKnIBXaget9jYWEaPHs24cePYsmULnTp1omfPnsTHxxda/5tvvqFbt27ExcWxefNmOnfuTO/evdmyZYu7jt1up1u3bhw4cID333+fXbt2MW/ePOrWrVvSj+O86eSEYqpXrx5r1nxFgwYN+N//3mbcuDHeDklEiujEiVPUrFnD22GIVCqukxOKl4kZmMV+zssvv0x0dDRDhw4FYMaMGaxcuZLZs2czbdq0AvVnzJiR7+upU6fy4Ycf8vHHH9OiRQsA5s+fz4kTJ9iwYYP7vPKoqKhix+ZJ6nErgVq1atGnTx+++OIbb4ciIsWQnHyS0NDQf64oImVCampqvis7O7vQena7nc2bN9O9e/d85d27d2fDhg1FepbT6SQtLY0aNc7+cffRRx/Rrl07RowYQXh4OM2bN2fq1KnuLYa8QYlbCZimyffff0d4uObJiJQnJ04ocRO50M5nqDQyMpKQkBD3VVjPGUBSUhIOh4Pw8PB85eHh4SQmJhYpzunTp5Oenk7fvn3dZfv37+f999/H4XAQFxfH008/zfTp05kyZUrJPgwP0FBpCcXHHyI6ur+3wxCRYjhx4hR16nh3mEOksinJYoMz9Q8dOkRwcLC73M/P72/vM/5yxqlpmgXKChMTE8OkSZP48MMPCQsLc5c7nU7CwsKYO3cuVquVVq1a8ccff/Diiy8yYcKEor8hD1LiVgKGYRARUZt9+w56OxQRKYbk5JNcfnlLb4chUqmUZLHBmcQtODg4X+J2LqGhoVit1gK9a8eOHSvQC/dXsbGxREdHs3TpUrp27ZrvexEREfj6+uY7Mq9JkyYkJiZit9ux2WxFe0MepKHSErrjjjt5//2PycoqfLxdRMoW0zQ5ciTBvR+jiFwYF2IDXpvNRqtWrVi9enW+8tWrV9O+fftz3hcTE8OQIUNYvHgxvXr1KvD9Dh06sHfvXpxOp7ts9+7dREREeCVpAyVuJXb//fdz8uQpPv644J4vIlL2HDuWxKlTKTRp0sTboYhUKhdqO5BHHnmEN954g/nz57Njxw7GjBlDfHw8w4cPB2Ds2LEMHjzYXT8mJobBgwczffp02rZtS2JiIomJiaSkpLjrPPjggyQnJzNq1Ch2797Np59+ytSpUxkxYsR5fy4lpcSthC677DI6dOjA8uXe24RPRIpu9+59ADRt2tTLkYhIaejXrx8zZsxg8uTJXHXVVXzzzTfExcW5t+9ISEjIt6fbnDlzyM3NdW+ue+YaNWqUu05kZCSrVq3ihx9+4IorrmDkyJGMGjWKJ5988oK/vzM0x+083HjjjUyf/hIOhyPf+LeIlD3ff7+ZgIAAGjRo4O1QRCoVI+8q7j0l8dBDD/HQQw8V+r0FCxbk+/rrr78uUpvt2rXju+++K2FEnqcet/Nwww03kJKSyi+/bPd2KCLyN5xOJ++88wF9+/Z1b6IpIheGzir1LCVu5+Hqq6+mSpUqfP110Tb3ExHvWLt2AwcPHmLYsGHeDkWk0rkQixMqEyVu58Fms9GrVy8WLnyP3Nxcb4cjIufw9tuxNG3alHbt2nk7FJFKRz1unqXE7TyNHTuWAwfiWbAg1tuhiEghjh1L4tNPv2DYsGFF2ohTRKQsU+J2nq666iruv/9+nnnmRQ4d+sPb4YjIX3z0kWvLnoEDB3o5EpHKyXXIfPEu/Yl1bkrcPODFF18kJKQaY8aMxzRNb4cjIn+yatVXdOrUiZo1a3o7FJFKySjhJYVT4uYBISEhzJkzhy+//IbY2BXeDkdE8uTk5LBu3Xf07NnT26GIVFoWwyjRJYVT4uYhvXr14u6772bs2CkcPXrc2+GICLB9+26ysrL/9sgbESld6nHzLCVuHjRjxgx8fW088cRkb4ciIsBPP/2C1WqlRYsW3g5FpNJS4uZZStw8KDQ0lNdee40PP/xcZ5iKlAFbt/5Ks2ZNCQwM9HYoIiIeocTNw+6880769OnDY489w4kTJ70djkil9tNP27j66mu8HYZIpaYeN89S4uZhhmEwa9YscnMdDBw4gszMLG+HJFIpZWRksmPHbq6++mpvhyJSqRmGUaJLCqfErRRERETw6aef8vPPv3HffaPIycnxdkgilc6vv+7E4XAocRPxMvW4eZYSt1LStm1bli1bxpdfrmPQoBEcP57s7ZBEKpUtW37BZrPRvHlzb4ciUqkVd/PdM5cUTp9NKbrxxhtZvnw5mzdvo337Xnz//U/eDkmkUjBNk08//YKrr26NzWbzdjgilZphlOySwilxK2W9evXit99+o2nTZtx22xC+/HKdt0MSqfC++mo969Z9x+OPP+HtUEREPEqJ2wUQFhbG559/zg033MBddz3A8uVx3g5JpMJyOp1MnjydDh3a07t3b2+HI1LpGSX8nxTOx9sBVBYBAQEsW7ac++67j+jo0aSmpnHPPf28HZZIhbNixWf8/PNvrFu3TivTRMqAkiw20E/uuSlxu4B8fX15++23qVatGqNHP82RIwk8/vj/4eOj/xtEPOHEiZNMmPA8N93Ui44dO3o7HBFBiZunKWO4wCwWCzNnzqROnTo8/fTTrFv3PW+88TJ160Z4OzSRcs1utzNixJNkZmYze/b/vB2OiOSxAJZiZmIWs1RCqRA0x80LDMNg7NixrF27lsOHExkw4EHsdru3wxIpt06cOMltt93LmjXrWbhwIRdddJG3QxKRPJrj5llK3LyoY8eOfPjhh+zYsZupU1/1djgi5dKePfvp1u1Odu3az5o1a+jVq5e3QxIRKTVK3LysZcuWPPvss8ycOY9vv/3e2+GIlBumafLhh5/Rvfud+PtX4fvvv6dDhw7eDktECqFTEzxHiVsZ8Nhjj9GpUycefPBxUlJSvR2OSJm2bdsOnn32ZTp3vpUhQ0Zy/fWd2bBhAw0aNPB2aCJSCG3A61lK3MoAq9XKwoULSUk5zb///Yy3wxEps9au3UD37nfyzjsf0LBhI1auXMny5SsICQnxdmgicg46q9SzlLiVEVFRUfz3v/9l6dKPNGQqUoivv/6W/v2H0blzZw4ePMgHH3xA9+7dvR2WiPwDC0aJLimcErcyZNCgQbRu3ZpJk17CNLUWWuSMr75az113PcANN9zAsmXL8ff393ZIIlJE6nHzLCVuZYhhGDz33HP8+ONWPv10tbfDESkT1qxZx113DVfSJiKCErcyp0uXLnTt2oX//OcVcnNzvR2OiFd9+eU6Bgx4kK5du7Bs2XL8/Py8HZKIFJMWJ3iWErcy6LnnnmfXrr0sWbLc26GIeM2aNeu4++4H6datKx98sExJm0g5paFSz1LiVga1atWKO++8k5demq25blIp7dt3gHvvHUWXLjfw/vsfKGkTKcd0coJnKXErox566CEOHjzEjz9u9XYoIhdUenoGgwePICKiDjExS5S0iZRzFqNklxROiVsZ1alTJyIiIli27FNvhyJyQf3735M4cOAwH3zwAcHBwd4OR0TO04UcKp01axb169fH39+fVq1asW7dunPWXbZsGd26daNWrVoEBwfTrl07Vq5cec76S5YswTAM+vTpU8LoPEOJWxlltVrp27cvK1Z8hsPh8HY4IhfEtm07iIlZzowZM2jWrJm3wxGRciQ2NpbRo0czbtw4tmzZQqdOnejZsyfx8fGF1v/mm2/o1q0bcXFxbN68mc6dO9O7d2+2bNlSoO7Bgwfdpxx5m2EWYRJVamoqISEhpKSk6C/gC2jjxo20b9+ejz9+h44d23g7HJFSFx09ms2bf2XPnj34+vp6OxyRcs+b/36fefbnoRFUsRSvnyjd6aRHUkKx4m7Tpg0tW7Zk9uzZ7rImTZrQp08fpk2bVqQ2mjVrRr9+/ZgwYYK7zOFwcN1113Hvvfeybt06Tp06xYoVK4r1fjxJPW5lWNu2bYmKiuKDDz7xdigipW7//oOsWPEZjz/+uJI2kQrkfBYnpKam5ruys7MLfYbdbmfz5s0FTlPp3r07GzZsKFKcTqeTtLQ0atSoka988uTJ1KpVi+jo6BK8e89T4laGGYZBv379+OijleTk5Hg7HJFSNXPmPEJDQ7n33nu9HYqIeND57OMWGRlJSEiI+zpXz1lSUhIOh4Pw8PB85eHh4SQmJhYpzunTp5Oenk7fvn3dZd9++y1vvvkm8+bNK9mbLwVK3Mq4/v37c+LESdau3ejtUERKzUcffc7bb8fy5JNPEhAQ4O1wRMSDLCW8AA4dOkRKSor7Gjt27N8+y/jLzr2maRYoK0xMTAyTJk0iNjaWsLAwANLS0hg4cCDz5rn+qCwrfLwdgPy9q666issuu5Rlyz6ha9drvR2OiMdt3vwzDzzwGP369WPUqFHeDkdEPKwkq0TP1A8ODi7SHLfQ0FCsVmuB3rVjx44V6IX7q9jYWKKjo1m6dCldu3Z1l+/bt48DBw7Qu3dvd5nT6QTAx8eHXbt20bBhwyK+I89Rj1sZZxgG/fvfxaeffkFGRqa3wxHxqPT0DIYMGclVV13FW2+9haWYE5hFRABsNhutWrVi9er853yvXr2a9u3bn/O+mJgYhgwZwuLFi+nVq1e+7zVu3Jht27axdetW93XzzTfTuXNntm7dSmRkZKm8l3+iHrdy4J577uE///kPMTHLiI6+29vhiHjMK6/8j+PHk1m7dp2GSEUqKsMo0nDlX+8prkceeYRBgwbRunVr2rVrx9y5c4mPj2f48OEAjB07liNHjrBw4ULAlbQNHjyYV199lbZt27p76wICAggJCcHf35/mzZvne0a1atUACpRfSPrzthxo0KABt99+O6+//pb2dJMK4/ffD/Lf/77Jv//9bxo0aODtcESklFyoDXj79evHjBkzmDx5MldddRXffPMNcXFxREVFAZCQkJBvT7c5c+aQm5vLiBEjiIiIcF9lfcqG9nErJ3744QeuueYa3n77v9x8cw9vhyNy3gYMGM6vv+5m586dBAYGejsckQqpLOzj9lVYXaoWcxrEaaeTzseOKO8ohHrcyomrr76a6667jpkz39DB81LuffHFN3z22ZdMnz5dSZtIBWfkDZUW95LCKXErRx5//HE2b/6ZL78899lrImVdWtppHntsEp07d+aOO+7wdjgiUsp0yLxnKXErR3r27EmnTp14+ulp2pBXyq3x458jKekEb7zxhv6qFhEpJiVu5YhhGLz66qvs3r2Pt96K8XY4IsW2evVa3n47lunTp2tBgkglYViMEl1SOCVu5UyLFi247777mDZtJidOnPR2OCJFdvLkKUaNGkePHjcybNgwb4cjIhfI+Rx5JQUpcSuHpkyZgsPh5Lnn/uvtUESKxDRNxowZT1aWnTfeeFNDpCKViBI3z1LiVg6Fh4czfvx45s9fzI4de7wdjsg/eued9/nww8+ZN28edevW9XY4InIBaVWpZylxK6dGjhxJ/fr1GTduirYHkTJtz579PPnkswwdOpTbb7/d2+GIyAWmHjfPUuJWTvn5+TF9+nS++upbPvtsjbfDESnUqVMpDB78f0RGXsyMGTO8HY6ISLmnxK0c6927Nz163Mhjj03UQgUpczIzs7jrruEcP36CFStWUKVKFW+HJCJeoKFSz1LiVo4ZhsEbb7xJdnYOI0eO05CplBm5ubkMHTqGn3/+jU8++YTGjRt7OyQR8RINlXqWErdyrm7durz55pt8+ulqFixY4u1wRDBNk8cem8TKlV+xdOlS2rZt6+2QRMSLLIZRoksKp8StAujTpw8PPPAA48ZNZdeuvd4ORyox0zT5z39e4e23Y3njjTfo1auXt0MSES9Tj5tnKXGrIF5++WXq1avP/fc/SmZmlrfDkUroxImTDBz4EC+/PJvnn3+eIUOGeDskESkDDEowxw1lbueixK2CCAwMJCYmhr17f2fUKM13kwsnJyeHV1+dyzXX9OC7735ixYoVPP74494OS0SkQlLiVoFceeWVLFiwgKVLP2LGjLneDkcqgeTkE9x227385z+vcMcdd7Jt2zZuueUWb4clImWIYSnZJYXz8XYA4ll9+/bl119/5dln/0OjRg3517+6ejskqaC2bdvOwIEPkZlp58svv+Taa6/1dkgiUhaVZHsPTXI7J+W0FdCkSZO47bbbGDbsUeLjj3g7HKmAPv10NT169Cc0NJwff/xRSZuInJMWJ3iWErcKyGKx8NZbb1G1ahAvvviat8ORCmbRoqUMHvx/9Op1E+vWrePiiy/2dkgiUoa5ErHiLlDwdtRllxK3CiooKIgHHniATz5ZjcPh8HY4UkHMnDmPkSOf4v777ycmJobAwEBvhyQiZZx63DxLiVsF1rVrV06dSuGXX7Z7OxQp50zTZOLEF5g48QXGjRvH7NmzsVqt3g5LRKTSUeJWgbVp04bAwEDWrt3g7VCkHHM4HIwe/TQzZ87jlVde4T//+Y/OERSRItPJCZ6lxK0Cs9lsXHfdtXzzzUZvhyLllGmajBkznnff/YAFCxYwevRob4ckIuWMhko9S4lbBdelS1c2bvyRrKxsb4ci5YxpmkyY8DyLFi1l/vz53HPPPd4OSUTKoeIvTCjB9iGViBK3Cq5NmzZkZWWzf/9Bb4ci5cz06bN57bU3mTlzJoMHD/Z2OCJSTqnHzbOUuFVw9evXB+DAgUNejkTKk7lzFzFlyis888wzPPzww94OR0TKMSVunqWTEyq4atWqAZCenu7dQKRccDgcTJr0Iq+99iZjxoxh/Pjx3g5JRET+RIlbBZeTkwOAzebr5UikrDt1KoXo6NF8/fUGXnnlFUaNGqV5JiJy3gyLgWEp3u8Sw9TvnnNR4lbBnUncfH2VuMm57dy5h7vvfoiTJ1NYuXIlXbvqjFsR8YySDH3qb8Zz0xy3Cm7VqlUAREVFejkSKas+++xLunfvS0BAFX744QclbSLiUdrHzbPU41aBpaamMnnyM9x4Y2eaNWvk7XCkjElNTWPatFeZM2cht9xyCwsXLiQoKMjbYYlIBVPZe9wOHTrEgQMHyMjIoFatWjRr1gw/P78St6fErQLZtm0b48aN45ZbbqFZs2YMHjyIo0eP8uabL3s7NClDTNNkxYrPGDduKikpabz44ouMGTMGi0Ud8CLieSXZl628z689ePAg//vf/4iJieHQoUOYpun+ns1mo1OnTgwbNozbb7+92L979Zu6gkhOTub6669n48YNDB06lHbt2mGz+bBmzTKaN2/s7fCkjDh9Op377hvFffeNok2bduzYsYNHH31USZuIVAizZs2ifv36+Pv706pVK9atW3fOusuWLaNbt27UqlWL4OBg2rVrx8qVK/PVmTdvHp06daJ69epUr16drl27smnTpr+NYdSoUVx++eXs2bOHyZMn89tvv5GSkoLdbicxMZG4uDg6duzI+PHjueKKK/jhhx+K9R7V41ZBrFixghMnTrB9+3p++ukXqlQJpEOHa7QoQdwOHIjn7rsfIj7+CO+99x533nmnt0MSkUrAoARDpSV4TmxsLKNHj2bWrFl06NCBOXPm0LNnT7Zv387FF19coP4333xDt27dmDp1KtWqVeOtt96id+/efP/997Ro0QKAr7/+mrvuuov27dvj7+/PCy+8QPfu3fntt9+oW7duoXHYbDb27dtHrVq1CnwvLCyMG264gRtuuIGJEycSFxfHwYMHufrqq4v8Pg3zz/1355CamkpISAgpKSkEBwcXuXG5cNasWUOXLl3YsCGOJk0u9XY4UsZ89dV67rtvNKGhtVixYgXNmjXzdkgicgF489/vM8/+vXVjgnysxbo3LddB/R93FivuNm3a0LJlS2bPnu0ua9KkCX369GHatGlFaqNZs2b069ePCRMmFPp9h8NB9erVee2117x2oozGRyqIq6++mmrVqvHYYxPJzMzydjhSRjidTl577U3uuCOatm3bsmnTJiVtInJhleTUhGJ2udntdjZv3kz37t3zlXfv3p0NGzYUqQ2n00laWho1atQ4Z52MjAxycnL+ts6fZWZmkpGR4f764MGDzJgxo8CQbHEocasggoKC+PTTT/nxx595660Yb4cjXmaaJu+99yFt2vRg/PjneOyxx/jkk0+pXr26t0MTkUrmfA6ZT01NzXdlZ2cX+oykpCQcDgfh4eH5ysPDw0lMTCxSnNOnTyc9PZ2+ffues86TTz5J3bp1i7xt0pkV+wCnTp2iTZs2TJ8+nT59+uTrGSwOJW4VSPv27enRowfz5y9m374D3g5HvCQh4Sh33fUADzzwGM2bX8mGDRt4/vnnsVqLN1QhIuJtkZGRhISEuK9/GvL862pU0zSLtEI1JiaGSZMmERsbS1hYWKF1XnjhBWJiYli2bBn+/v5Fiv+nn36iU6dOALz//vuEh4dz8OBBFi5cyMyZM4vUxl9pcUIF8+yzz3LbbbfRrdudbN68murVq3k7JLlATNNkyZLlPPXUFPz9A1mxYgW33HKLt8MSkUrOsLiu4t4Drj3Q/jzH7Vz7n4WGhmK1Wgv0rh07dqxAL9xfxcbGEh0dzdKlS8/Zk/bSSy8xdepUvvjiC6644ooiv4+MjAz3/pirVq3itttuw2Kx0LZtWw4ePFjkdv5MPW4VzBVXXMG3336L3Z7DK6/M8XY4coH88Uci/frdz0MPPcFNN93Mb7/9pqRNRMqE8xkqDQ4OznedK3Gz2Wy0atWK1atX5ytfvXo17du3P2dsMTExDBkyhMWLF9OrV69C67z44os8++yzfP7557Ru3bpY7/2SSy5hxYoVHDp0iJUrV7rn4B07dqzEi0XU41YBhYeH8+ijjzJ16lR8fHz4979HEBBQtG5dKV8yM7OYPXsBM2bMoUqVqnz44YfcfPPN3g5LROQsi+G6intPMT3yyCMMGjSI1q1b065dO+bOnUt8fDzDhw8HYOzYsRw5csQ95ywmJobBgwfz6quv0rZtW3dvXUBAACEhIYBreHT8+PEsXryYevXquetUrVqVqlWr/mNMEyZMYMCAAYwZM4YuXbrQrl07wNX7dmbLkeJS4lZBjRs3DovFwtSpU8nKymLq1HHeDkk8yDRNPvjgEyZPfomEhGM8+OCDPPPMM1p8ICJlzwU686pfv34kJyczefJkEhISaN68OXFxcURFRQGQkJBAfHy8u/6cOXPIzc1lxIgRjBgxwl1+zz33sGDBAsC1oa/dbueOO+7I96yJEycyadKkf4zpjjvuoGPHjiQkJHDllVe6y7t06cKtt95a7PcI2setwnv22WeZMmUK330XR716BTcglPIlNTWN99//mLfeWsKvv+6gT58+PP/881x22WXeDk1EyqCysI/boU6XE1zMfdxScx1ErttW7vOOBQsW0K9fPwICAjzWpua4VXCPPPIIYWFhPPHEsxQhR5cyasuWbYwaNY6mTTvy738/Q8OGl/H111+zfPlyJW0iImXU2LFjCQ8PJzo6usj7yf0TDZVWcFWqVOHVV1/ltttu47PPvuRf/yra3jNSehwOB0eOJLB//0GOH0+mdu0w6tSpTUREOIGBZ/8qO3ToD1au/Ip33lnKzz//RmRkJI8//gTR0dHnPGpFRKTMuUBz3Mqiw4cP8+mnn7JgwQI6d+5M/fr1uffee7nnnnuoXbt2idrUUGklYJomPXrcyNGjf/D11yu8HU6l9tNPv3D//Y+yf/+BQr9fvXo16tSpTU5ODrt378NqtdKzZw8eeGA4PXv21F5sIlIsZWKo9PorSzZU+vXPFSrvOHbsGO+88w4LFixg586d9OjRg+joaHr37o3FUvQBUPW4VQKGYTBgwN0MGTKElJQ0QkKCvB1SpfT227E89tgkWrZswWuvzeKyyy4jLCyMxMREDh8+nO9yOBxMmfIcXbt2pVq1at4OXUSkxAyLgVHMHrTi1i8PwsLC6NChA7t27WL37t1s27aNIUOGuA+4v/7664vUjhK3SuDAgQMsX74cgBMnTihxK0WmaeJ0OgvtGXv++f/Sp08f3n33XWw2m7s8KCiISy+99EKGKSJy4VygVaVl1dGjR1m0aBFvvfUW+/fvp0+fPnzyySd07dqVzMxMnn76ae65554ib8irxK0SeOihB/nxxx957rmniYqK9HY4Hrdr1162bNnmPtqkWbPGNG/e+G+POdm+fTdLl37E3r37qVYthIsvrkuDBvVo0CCKBg2iCAkpete80+lkypQZLFmynOTkk+Tm5tK8eROuvvoqWre+ipYtL6dGjWocO5ZEly5d8iVtIiIVnWGUoMetgiRuvXv3ZuXKlVx22WXcf//9DB48ON8B9QEBATz66KO88sorRW5TiVslsG3bNgYNupNBg/qydOnHHD78B+npGTRqdAm33NIDf//Cd6Iuq778ch0rVsRx9Ohx9u8/WOi5rLVrh9G2bStCQoIJCqpK1apVqFq1CklJJ1i9ei2//baTGjVq0Lp1a3btOsDnn39NUlKS+/6aNWvQoEEUTZtexvXXd+C669oVenyYaZo89NATvPfeh/zf//0fl156KYZh8OOPP7Ju3UbeeOOdfPWbNGni6Y9DRETKqLCwMNauXeveeLcwERER/P7770VuU4lbJXD6dDoBAf7ccUc0Gzf+QM2aNalSJZD4+EOsWBFHTEz5OBorJSWVceOm8u67H9CsWTMaNmxIjx69uPHGG+ncuTN+fn7k5OSwceNGPv/8czZv3syhQ0dJS0vLu05TvXo12rVrz9Spz9OjR498vV+nTp1i37597Nmzh71797Jnzx42bdrE22/HYhgGLVteQefOHbj++vYEBwdht+fw1VffEhu7gsWLF3PXXXcViPnEiRP8/PPPJCUl0bJlSxo2bHghPzIREe+rxEOlb7755j/WMQzDvUlwUWhVaSXQrl1b9u3by/HjycTFxdGzZ0/Adf7ahAnjOXz45zK/WnHlyq8YM2Y86emZvPLKK9x7770XrCv90KFDrF69mlWrVrF69WpOnDiR7/sjRozgtddeuyCxiIgUR1lYVXrkxlYE+xavnyg1J5e6KzdXiLwjPT2dtWvXEh8fj91uz/e9kSNHFrs99bhVAvfcM4QHH3wQgMaNG7vLIyMjycrKJj09g+Dgsrlg4dSpFJ56agoxMcvp0eNG5s17g4suuuiCxhAZGcl9993Hfffdh8Ph4LfffsNut+Pn54e/vz+XXHLJBY1HRKQ8+fOh8cW5pyLYsmUL//rXv8jIyCA9PZ0aNWqQlJREYGAgYWFhStykcNHR0e7E7c97xeTk5ABw882DSE4+SZUqgfTv34fRox8olTi2bv2Vb7/dxP79B/HxsdKy5ZVcc81VhIXV4tSpFHbs2MPOnXtITU0DwOFwEhOzjIyMLObPn8+QIUO8/sNstVq54oorvBqDiEi5Uok34B0zZgy9e/dm9uzZVKtWje+++w5fX18GDhzIqFGjStSmErdKwNfXl+PHjxMTE0Nk5NlVpR07dqR7924EBQVz002N2bt3L8888xIvvPAaHTu25c03XyEoqGqJnnn06HE2bPiBY8eOc/RoEj/99Atr124gMDCQSy5pSFZWFnPnLipwX5UqVahZ8+yKm7Zt2zFz5n8veC+biIh4SCWe47Z161bmzJmD1WrFarWSnZ1NgwYNeOGFF7jnnnu47bbbit2mErdKIjQ0lIcffjhfWf369Vm5cpX7a9M0uemmmzh48CAvvPACTz75LK+//nyxnnP6dDpTp85g/vwYsrOz8fPzIyKiNlFR9Xjvvfe47bbb3PPpTpw4waZNmzh58iRVq1alefPmREVFFWsHaRERkbLK19fXPVIUHh5OfHw8TZo0ISQkhPj4+BK1qcRN3AzDYODAgQDExsayePEyxowZziWX1C/S/d9++z0jRozl2LEkxo0bxwMPPECtWrXOObxZo0YNevTo4bH4RUSk7DEsrqu491QELVq04Mcff+Syyy6jc+fOTJgwgaSkJBYtWsTll19eojYryEcjnrZv3z4A2rX7F5MmvUh6esY566anZ/DEE5O56aaBREZG8csvvzB+/HjCwsK8PidNRES87MxQaXGvCmDq1KlEREQA8Oyzz1KzZk0efPBBjh07xty5c0vUpnrcpFCJiYnY7XZmzZrF1KlT+eCDT3jwwSEEBVXFarVQu3YYV17ZjJ079/Lww0+RmHiMV155hZEjR2qoU0RE3CrzWaWtW7d2v65VqxZxcXHn3aYSNylUUJBre5Dx48dz9913M3r0KCZMeB6Hw1GgbocO7Vm5crXO2xQRkYIq8eKE0qDETf5RgwYN+OijjwHXAobc3Fz27dvHb7/9hmma3HrrrWVuA18zIR4zPRWjQRMMiys2055N5o8/4lv3Inz/tEu18+gfOI/EY216FYa/v6uuacLxw+Bjw6gRfrbdzNNw9BCER2IEnF1xa2akgtMBVaq5h4edGRlk/rQZW/0G+Nat666btGcfp+IPUa9jO3z8XMeNOXNzObLhOwJCQwltenavPWfScXL37san+RVYqpbNvfZERP5WJdsOpEWLFkWeJvTTTz8Vu30lblIshmHg6+tL48aN823mW5Y4N6/FuXgmAMYV7bDe8xim08mRewaS/fNWsFiImDufwPYdcOzcRuYTwyE3F0tUQwJefRvD1xfn95/Bvp9dbbTujqVRK8y0UzgXPQdZGRBQBcugJzGqhmAeO4h5aKfr4aEXYUQ1w2nP5lC/27Hv2QO+vkS+uwT/K65kz+o1LOh9J6bDwcXt2zDsq8+wWCx8PCiavR9+AsC/3ppDk3534Ig/wKkhfSEzA0t4bUIWLcMSVL53EBcRqej69Onjfp2VlcWsWbNo2rSp+7zS7777jt9++42HHnqoRO0rcZMKx/njWvdr85eNmLk55CYkuJK2PKc/jyOwfQdy168Bp9N138F9OOP3Y23YCH7/9Wwb+3+BRq0w43e5kjaAzHTM+F0YTa/BTDpy9uHJf0BUM+y7druSNgCHg9OrVuJ/xZX8EvsB5J0yF7/he1IOHSYoPMydtAHsWLKUJv3uwL7+a8h0Pc95NJHcX7Zg63CdBz8pEZHSV9lOTpg4caL79dChQxk5ciTPPvtsgTqHDh0qUfuaRS4VjqVB07Nf1K2P4eOLT3g41vBw17wJpxP/Fi0BsDa90pW4GQYEh2CJyBvSrPWnDX/DLgbACL/47Bp1i8X1NUDQ2Q2DqVoNAN969bBUq/an57UAIKpDW8y854VE1iUoojZWPz9Cmzdzz+mo264NAD7Nrzw7z8PPH+sll3nmAxIRuZDODJUW96oAli5dyuDBgwuUDxw4kA8++KBEbarHTSoco8ttWGqGw+kUjNauHirD5sdFS97n9Kef4FuvPlVu6AKAT/vr8Z/8Ks4De/Dp0AUj0DVvzXLd7Zj7t7nmuNVv7mojNALLXY+4etqiGmPUrO0qv+gyCAgCZy7UdCV+1qAgLn5/OadXfo6tcROqdOgIQKshAwmsWYPkffu5st8d+NhsAPT9/EO2L44lMKwWje907aTte0ULgme/Te4vW/DteD3W8IgL9AmKiHhSSbb3qBiJW0BAAOvXry+weG/9+vX4582pLi7DNPPGbf5GamoqISEhpKSkEBysOTYiIiLlgTf//T7z7GP9riXYVrx+olR7LmGx35T7vOO5555j0qRJDB06lLZt2wKuOW7z589nwoQJPPnkk8VuUz1uUuGYjlzMxP2Qa8cIi8IIcK3GNA/sxLnpC4wa4RjX98Hw8cV0OHCuWYaZcBBLm65YGl3lqpt4CMcn74JfANZbh2BUDQEg/YNYstavxb/T9VS5ra+rbnoa5tfLIScbo1Nvd0+ceWAHzu2bMGrVwWh1A4bFipljx/n5EsykBKzX3oTRsBkAubu2k7nwDSw1QgkYPhJLlaqYpon51Yc492zD0upaLK01v01EyqFKtqr0z5588kkaNGjAq6++yuLFiwFo0qQJCxYsoG/fviVqU4mbVDhmwj448YfrdfopaNIB7Fk4Y2ZAbi4mgJ8/RqfemN9/gblmGQDOXVswnngNI7g6ua9PguOJADgy0/EZ9hTZP3xPytRJAGR/8xU+9erj1/JqnJ+/A7u25D37INYH/4N5OgXnZwvBdGIe2A4BVTGat8NcswLz6w/BBMfOLVifeQtsfqSOfgAzJQUMMJ0Oqv57PObPG3Auf9MV26+bMOpEYdSpd8E+RxERKZndu3dz2WWuecl9+/YtcZJWGC1OkIonJ/vs69wcwITszLOvDQPSTgFgpp08u+DA6YSM067XKSfAdIJpYp5MAsCRdDzfY5zHj7lepJ50rRQ1TTid4irLynDdD67npae6npeSjGvuhgk59ry4cjFTU9zPc7ebciLf88zUk+f3uYiIeMGZVaXFvcqzFi1a0KRJE5544gk2btzo0baVuEmFY4RFgcXVmWzUboBhWDCCa8A1rgUJVAnCuKYbAJY2XaFaTVfdltdCuGs1qeX2oWCxgJ8f1psHAeB//Q34XnEVAL5XtsDv2htcda/rA74210rTLne6nlGzNlzmWrlK1WoYTdvk1e0NVV3zNYxrb8IIqYFhsxE47GEwDIyqQQTcc7/r+62vhzp5GwU3bYVxackOJBYR8apKuKo0OTmZF154geTkZG699VbCw8OJjo7mo48+Iisr67za1uIEqZDMvN6rM6cmuMvt2eDji/Gn81RNpxNy7Bh+/gXrWq0Y1rMzCkzTxMzIwAgMzPcXoenIdT3Pxzd/Gzl5zzP+/DwH5OQUfF5WpquuT/7nkZ2J4R9Ygk9BRCq7srA44figG0q0OKHWojUVIu8wTZONGzfy0Ucf8dFHH3Hw4EG6du3KLbfcwk033URYWFix2lOPm1RM9mw4nVawPDsTHLn5y3JzcaakUOBvGHu2azjzz0wnRvafhkHd5WaBsjNJ15kNft3lObk4UgvGZmZlFnye0wlZGa7kUkSkHDpzyHxxr4rCMAzat2/Pc889x/bt29m6dSvXXnstCxYsIDIyktdff71Y7WlxglQ4zu1byHl5HNizsd48EJ877sU0TZzvzcLcvBYCqmAdPgkjIorc+IMcGzIA54lk/Dt3oeb0/2JYLDhWvY/z/TfAxwfrA09jubItZkY62RNHYMbvx7i4IX6TX8cICMQ8fhjnb+tdSV3DFlguboJpOnGumAd7f4EqwVjuegSjei2yd+7g8OC7caamEnz7HYRPeQ6AzHkzsb/7Jvj5U+X51/G5qjVm6kkcL/8bkhKhQVOsD0/B8PX9h3cvIlLG6JD5fC699FIeffRRHn30UZKTkzlx4sQ/3/Qn6nGTCif3s6XunivHx+9i5ubAqWRX0gaQlYlzw0oA0pe/j/OUa9J/1ldfkrPXdUyV86NFeY3l4ohb4mpr87eY8fsBMOP34dj8ravuwV/dvW3mgW2u+44dcSVtABlpmNs2AHDq3UU4T7sWQKR+8D45iQmYdjv2xfNdde3ZZMe+7Wpr8zeupA1g/3bMPb947DMSERHvOXnyJAsXLqRmzZoFNuf9J0rcpMIxauWdMGCxuBYeWH0goArY/FwrSE0nRs1wAKx167qGIy0WsNmwhoa67q0Z7iozLBhhdfK3m/eXoBGW93VAEK6Vogb4V3GVVQkGq9VV1zShmqtd34siXV9bLBhVqmANDgEfH4waoa7nAZY6kXkx1D77PMPAqFG8eRAiImWChRIsTvB20KUrPj6ee++9t0T3aqhUKhyffkNxBARipp7E2rOvaxGBfwDWYRNxblyJUSsCo1MvAKrc1hczM5OcXTupcusdWGu4Vpj6PDwZ52ex4B+IpdcAAKyNL8d35AScW77D0qId1stcR2FZLrsa0y8AHLkYUa4NdY2qIVju+D/MX7+DsIswLm8HQPX7hoJpYj94kGoDBmIJdC06qPLyXLJjF2KpURO/gUNdbVx+DZYBD2Pu+RWjZUeM2pEX7kMUEfGQynbIPLgWZvydtLRC5mAXkVaVioiIVFBlYVVp8v3dCbYVb35uqj2HmvNWFTvuWbNm8eKLL5KQkECzZs2YMWMGnTp1KrTusmXLmD17Nlu3biU7O5tmzZoxadIkbrzxxnz1PvjgA8aPH8++ffto2LAhU6ZM4dZbb/3bOCwWy98mn6ZpYhgGDoejyO/tDPW4SYVjOhw4130GaaewXPsvjJAaADgTj+D45nOM2hdh7dTd/UOV9lkcWTt2EHzzLfhdcgkAjlOnSF2yGKNKFUL69cew+bna3r0N546fsDRrhXFJ87zn5WLu+MF1xFaTqzH8XL1ojgP7yPniM6yXNsb3uq6uuqbJqeXLsR84QPU77sB28cWuuknHyVy+FEvNmgTccgeG1bWNSe5P3+Hc/jPWNtdivbTJBfoERUQ86AItToiNjWX06NHMmjWLDh06MGfOHHr27Mn27du5OO937Z998803dOvWjalTp1KtWjXeeustevfuzffff0+LFi0A2LhxI/369ePZZ5/l1ltvZfny5fTt25f169fTpk2bc8YSFBTEuHHjzllnz549PPDAA8V+j6AeN6mAcpcvwPlpjGs+W+2L8J08F3LsZI24E9JSwXTiO/RRfLrdQuonH5MwZjRYLFgCA2nw9TdYQ0I4PKAv2b/8DKZJcP8B1Bo/CfPQPhzPjXI/xzr2vxgX1cfx7Uew80fXL5qwSKw3DcWZmkJa/56Q5doOJGDiC9g630jyokX8MX4CWCxYa1Sn8fpvMfxsJN95E45DB8HppEr0g1Qd/jCOX38ia7xrY16sPgS8HoPlzLw6EZEiKBM9bg/0KFmP25zPixV3mzZtaNmyJbNnz3aXNWnShD59+jBt2rQitdGsWTP69evHhAkTAOjXrx+pqal89tln7jo9evSgevXqxMTEnLOdzp0707NnTx5//PFCv//zzz/TokULnCXY6qmCT/+Tysg8sDvvhRMS4sGRi3kqGVJPucosFpz7dwGQ9euvrkUBTifO06fJOXwYAPv27a5FC6ZJ9i+u1Zzmkd/PHm1lmq6vAY4fyXueCUkJrpcJRyAj3b3wwbF7OwCZv2xzP8+RlExuchLY7TgO/u7e7y1n528AOPftPttubg7OQ7+X5scmIlLmpKam5ruys7MLrWe329m8eTPdu3fPV969e3c2bNhQpGc5nU7S0tKoUaOGu2zjxo0F2rzxxhv/sc0BAwbg7+9/zu/Xrl2biRMnFimuv1LiJhWO9dqeZ1d+duiO4eOLUSsCy+WtXBUsVqzXuuYwBPe+GcNmA8D/yivxyzsUOLj/Xe72gu9yLU4wmrWGvGFXqtXEaOpqz2hyzdmHN7na9YiGl2Jt5FqogM2G7YaeAFS/43b3MGjVTp3wrVMHw88P/3/dnBebhYCbb3e9j7bXQtUg1zMiLsLa9Mrz/3BERC444+xwaVEvXL/DIyMjCQkJcV/n6jlLSkrC4XAQHh6erzw8PJzExMQiRTl9+nTS09PzHQifmJhYojbvv/9+Ro4cec7vh4eHlzhx0xw3qXAsrTrhO20BZnoaxsWuOWuGYWAb+yLm77sxaoS5tt8A/Js1o8HXa8k5dBj/pk3dG9zWfOIpgvrchiUgEN8o13mhRlA1rBPnwB8HoU6U+xgqS6NWmBH1XIfYV3f9gBs+vlR57W0ce3ZgiaiLpbprtWqVNm1otH4dOYmJBDRv7p5nFzxxKoF3DcYSEoI1oq6r3fA6BM5+D+fheCz1LylwRJaISLlgsbi3OyrWPcChQ4fyDZX6+fn97W1/XRBwZhHAP4mJiWHSpEl8+OGHBY6gKmmbpUWJm1RIRmhtjNDa+cusPhiXNC1Q16dmKD41Q/PXNQz8GhdcDGD4B0KDQsqDaxYs8/XFp+kVBcp9w8Px/ctfcIbFgm/jgrEZVYOxNm5eoFxEpNw4j8UJwcHBRZrjFhoaitVqLdATduzYsQI9Zn8VGxtLdHQ0S5cupWvXrvm+V7t27WK3uWTJEvr37/+PMYMrMY2Pj6dDhw5Fqg8aKpUKKHvfPvZ278au1q04tXy5uzxt4Xz+6HQ1x+6+A8exowA4T50kKXogCde2JvV//3XXdWz6hsz7byFrZH/3fDhME9KSIfmI679563rMU8dwrFqI47M3MRP2u9twblqFY94EHCv+h5mV4SpLPkbmI/eSflc3cj48O7E16b2lbGl6Jb9e14WsvftcdTMzOfLAMPa0uJKjkybovFIRKZ+KO0xagkTPZrPRqlUrVq9ena989erVtG/f/pz3xcTEMGTIEBYvXkyvXr0KfL9du3YF2ly1atXftjl79mwaN27M888/z44dOwp8PyUlhbi4OAYMGECrVq105JXIsekvYf/9dxwnT5Iw9knM3FwcScdJfeUFzNNp5OzcTtrbbwJwOvZd7L9swUxP5/S8WeQedC0AsP/veUg9iXksgZxFeQcA2zNdF2be6ywAnNvWQ0YqZGfi3LIGAPPkMcwfv4ScbEg44D7yKuf9RTj374b009jfnIl56gTOnBwOPvYkjpMnydq7jyPPvwRA6ocrSP9qDWZ6OimLF5P5w6YL+CmKiHjIBUjcAB555BHeeOMN5s+fz44dOxgzZgzx8fEMHz4cgLFjxzJ48GB3/ZiYGAYPHsz06dNp27YtiYmJJCYmkpKS4q4zatQoVq1axfPPP8/OnTt5/vnn+eKLLxg9evQ541i7di0vvfQSa9asoXnz5gQHB3PppZdy+eWXc9FFF1GzZk2io6OpV68ev/76K7179y7W+9RQqVQ4hn9A3gvDtfDAMMDHByxWcLo2OzTy5kkYfn7unjMMA/IWKmDzg4zTeWX+Z7+f70F5/7X6nD3ayuJztuwME/BxzZ07sx+cq47VddyVYWD4+mLm5gJgCXA9z/KXOW3u9yUiIgX069eP5ORkJk+eTEJCAs2bNycuLo6ovHnKCQkJxMfHu+vPmTOH3NxcRowYwYgRI9zl99xzDwsWLACgffv2LFmyhKeffprx48fTsGFDYmNj/3YPN4CbbrqJm266ieTkZNavX8+BAwfIzMwkNDSUFi1a0KJFCyzFnfeXR/u4SYWTc/QoCeOeIvf4ccIef5yqHToCkLEyjrT5c/Ct14BqT0/GEhSEMzOTlOcmk7NzO1UHDCbwFteKTseubeQsfA0joAq+wx7DElbHlZhlprp622wBEBAMhoGZnoJz61eQm4Pl8k4YNVxz65zbN2Fu+xZq1cVy7a0YPr6Y6afJfu05zD/i8b3zHnw6dgEg5etvODLteXzDwoh66Tls4eGYubkcf24aGRs3EHzrbdQYer93PlARKbfKxD5uI/sQ7FfMfdyyc6g5c4XyjkIocRMREamgykTiNurWkiVury5X3lEIzXGTCunUd5s49kkczj9t1piTdpr45R9yYusv+eo6D+3H8d1XmOl/OvQ3b9NbHLlnh1IBM/M0zoT9mJmn/1TVxLlzK86fv8P807lz5ulUcjd8ifPQ2QULAGbSHzj3/4ppPxubMyuLpE/jSNn0Q766jkMHyV71Kc7k4yX7IEREvO0CzXGrLDTHTSqcw2+9zc7HnwKgxnXX0vL9GJy5uazu/i9SduwEw6DTOwu4qFdPnDu2kvPC464TFWrVxjblTdc8NHuGK2kD8PUDX3/MjDScGz90ncRg9cXS/haMgKo4P4/F+dFCAIyrr8Pnvicws7PIfup+zOOJYFiwPf0K1qZXYR7YjvPLJQCYNcKx3PIghsXCbwOHcGrdegAuefE5IgYPJHfPLk7d1x9y7Bgh1ai+5GMs1WsUfMMiImXZBTqrtLJQj5tUOMfiVrpfn1j7DU67nYzDh11JG4BhcOTzVQA4f/7+7I3HEzH/yJu4eiZpA1fPG2CeTDxb7shxfQ04f1rvrmpu3ej67+EDrqQtj3NLXnn8rrO/kE4chfRUHJmZ7qQNIDkvNvumjZBjd92Xcorc7duK/VmIiHjdmQ14i3tJofTJSIUT2qWz+3W1dm2w2GwE1q1LUMMGrkKnk4i8OpZmLV29bQDVQzEiLnK9tvypM9qatyK0erhrZSqAxYpRzbW7tuXys0deuY/BqnsxVM/b1Nd0Yrm8tev1RZecHXoNCYUqQVj8/Qm6urW7jeqdrwPA1uoa18pTwKgahE8hG/SKiEjlosUJUuGYpknymq+xJyURflMvrFVcR1NlnzzJ4U8+I6hhA8Lat3XXd+7bgXn4AJar2mKEVD/TCDhyAOPsdh+AmZ6CefIoRvVwjCoh7ueZv3wH2VkYLTq6j80yTyXj2PIdlsgGWC45e9qCmXAAM+0ERlQTDD/XFh+O9AySPo3DViuUatdf5z5OJXfvbnJ//RnfNh2wRtQp1c9NRCqeMrE44bG+JVuc8NJ75TbvaNq0KevXr3cfWD9s2DCmTJlCrVq1ANfpC/Xq1SMjI6PYbStxExERqaDKROL2774E+9mKd2+2nZovlt/EzWKxkJiY6D73NDg4mK1bt9KggWvk5+jRo0REROAswYk4GiqVCsd58gTpE/9N2sj7yP3lJ3e548f1ZE0Ygf1/z7mPoDJz7GS/8TJZ44aTu/4Ld10zMw3n77/gjP8NM8e1+tM0TcyEfTh3bMRM2MeZv3mcxxLJmPgo6U88hGPvTncbW2KW8nqn7nw4+gly81a35pxO54eRY1jT82b+yJvLBpC8YSMbbr6dzfc/SHZysut5TieHnnuBX2+6hcS33i6lT0tEpJRpVSmF9ZGV9KB6rSqVCidz9svkfL0KTJP0J/6P4E/XQ0Y69ulPg8OBY9c2jGqh+PYfSu5nH+BYuRxME/vu37A0vhxLaDhm/HbXcVW4fuCMqOaQmoR5ZLerLP0URmAIhISSNXMauT9sABMyJj1G0DufcPJgPEvuGYbpNDn43Q9Ur3cx144ewY6XZ/D7u0vANPl28H3csmc7vlWrsOnuIeSmpYFhYA0I4KqZL5O84kOOzHCdn3r6x5+o2rIFVa8seGi9iEhZZlgsGMVcbFDc+pWJEjepcMy01LwXJmZmJjidriTszB5rhoGZ4dqzzUw/ffa4KtOErExXnT+vKnWvJP1TGeTNgQPzdJr7fjPdtb9b9ul0TKfrLyzDYpB1ynX2XU5qGoZhYDqdmDk5OLOzMAP8cWRkuNowDHLyzsnLTUnN/7i0NEREyp+S9KCV7x43wzAK9KiVtIftr5TSSoUTED0Co1Y4+PkRMGYsho8PRvVQfPreB1YfjIhIfHrfBYBvz9sxoi4BqxWf3v0x6rrOtDPqXOpaju7ji1E7bzVqtTAIqQUYEBIG1cIB8Lt/lGtRQ0AgASPHAhDetDEdRz2ExdeHOldeTvuHXMdVNR75EFUb1Mdis9F83BP4h4W5Xk/7DxZ/fwLqRNDoiccAqHXn7QS3bwtWK6F33EZwu7aIiEjZZ5omXbp0oWXLlrRs2ZLMzEx69+7t/rpbt24lbluLE0RERCqosrA44cTYuwn2L+bihCw7Naa9W27zjmeeeaZI9SZOnFjstjVUKhVS7vHjONLSsNWv7+6eNp1O7Pv24VOrFtZq1dx17SmpZBw5QnCjy7Dk7ZsGkH3gIBZ/f3xrh7vLzBw7HD0M4Rdh+J79RXT6jwQcWVmENKjvLnM6HCTs2EWNyLoEhIScbSM9DdJOQnjk2dhME/vvv2MNCsInb7k4gDMjg9zDh/CtVx/DVrxffCIiZUIlPDmhJAlZUWmoVCqc1C+/ZGe7duzp0oWEvB8e0zQ5/OBwfu/Zg70dO5C5dSsAp7bv4OPLW7Kyw/V83ecOnHnz4P544SW2d7iWX1u34cSKD11tnE4ld8L95E5+iNyJw9xnm+6KXcrbja7gnSuu5vv/PAeAIzeXV27oxbOXt+WpqGb8sd212tQ8sAvHxPtwTHsYx5vT3CuNEsaPZ2/Xruzq0IG0r74CIDchgUM9unL41t4c6XcbzhLs9yMi4nU6OaGAn3/+GeufOgqKo2J/MlIpnVi0yL0Q4cSiRTjtdnKOHCH9yy8BMO12Ti19D4CDS5biyEgH4Pj6DaRs3wHAsTlzXY2ZJsfnvel6uW0TJB91lSclur4Gfn7tf+7TELbOfB2AI7/8yp5vNgCQdfo03y+MAcC5YRXk5i1y2PY9nEzCmZ3NycWLXWUOByfeeQeA06s+x5GcBIB9926yftzkwU9JROQC0XYghSrCTLVCKXGTCse/USPXD73Vii0qCsPXF5+aNbGEhLiOrHI68bvsMgCCmzbGdDgxrFasVQIJrOs6ncD/sstcf/EZBgFN8049iLjY9V/D9WNj5H1ds3kzsFgwrFZqNGkMQI2LI7EFBmJYrZgOJ3Wau9ow6kS5jtgyLBAYBFWDMWw2fCMj3cdb+TdqBIDtUleMWCxgteIbdXYYVkSk3FDiVijt4yaSJ+zRR/EJDSU3OZkagwa5lmUHBBC1JJZTS5diq1ePav37A1Cv353gdHLq19+o178vfnnHkzR8ez7H31yApWoVwoZGA2Cpdxk8/Czm9s0YzVtjRF0KwLXTnyOkQX1y0tO58qEHAKgaWpPH1n3Od4uWEHnV5Vxzdz8AjGtvwuLji3k8AUu7bhg2P1ccixdzYtEifEJDqTF4MACB7TsQPuO/ZP74I1W6dcc3KurCfYgiIlImaVWpiIhIBVUmVpVOvLdkq0qfeavc5h2pqal/+/1ffvmF6667DseZ/UWLQT1uUuE4s7JImDOPnONJ1L7/PvzzeqrSftvOHwvfJaBBfSKH3usaxjRNkt5ZTOb27dTseydVWlwFgD0hkcT/zcVatQq1HxqOtUoVAJI/W8nJNV9TvUtnavboDkDu6XT2vvY6uenpNHzoQQIiagOQ9fNWUlcsw69xU4L79svbeNeBY+UyzMQjWLvcjOVi1x5xzsQj5MYtxQipgc/N/c+uWM21uzb69bGBtXiHNIuIlAklWWxQzhcnVKtW7W+HQk3T1FCpyBmHpr1A4tw3wGLh5MpVXLVpA86MTDbffBuO9AzXwgXTycXDh5Ecu5T4J54Cq4Wk2Pe54oeN+NSswe5BQ8jcuct1FNaRBOrPeIm0zVvYOWQoWCwcXfgOV3z+MUEtruKXJ8Zy6L33MQyDpPUbuP6r1eSeSObw4Lsxc3LA6cSw2Qi+9TYcn71P7ruzwWLB8e0X+L3+Ptj8yJ40EvPE8bzTF9KwDR4BuTmQ7Vo4Qa4dAoJdc/RERMqTSrgdyFd5uwOUBiVuUuFk7tnreuFwYD98BDMnB/uJEzhS846MslpJz6uTtXef6y87hxPTkYU9MRGfmjXI3rffvTI1c7frfNLM/ftd9zudrq/3/U5Qi6tI27kLnE5M4PTefQDkJiZi5h0sj8WCfb+r3Dxy0PU8pxMyTsPpVAgKwUzKW61qGDgPH3S9Nv/She50KnETkfKnEiZuBw8epF+/fvj5+Xm87fLdFylSiNpD78Xwcf1NUnvYUCw2G/4X1SWsdy8ArP7+1B08EICa/e7AGhQEQNC1nQho7FrRWfvhh1yNWa3UzltwUKNbF/zr1wPAv2EDanS7AYBLHn4II29F6KWjRwLg16gxAe07AGAJCiL41ttdzXW5GfKGQS1tO0ONWhg2P6z/uiPveT74ul/bzv7ysljBqr+zRKQcqoT7uN17772k5J077WlanCAVUu7JkzhOp+MXeZG7zDRNMg8cxBZaE5+8ZA3AkZFBTuJR/OpFYfzpl0X2kT+w+PnhG1rTXea028k+dBi/yIuw/Okkg+ykJJzZdgLythMB10kNOfHx+ISFYQkMPFuekY6ZdgojrE6+OQ7OYwkYAVUwgv70M2aaZ7cPKed/gYrIhVcmFidMeaBkixPGzSm3eYfFYiExMZGwsDCPt60/4aVCslarhjUk/w+7YRgEXBxZ4C85a2AglqiL8yVtALaI2gWSJYvNhl+9qHxHYwHYatYssJmiYbHgGxnp7o1zCwjE8PMvMDHVCA0vmJwZBqZplHgSq4iI1xmUYKi0VCK5oErr97YSN6lwnAf2kPvKOEhPw9r/Aaxd+wBweNoLJL42C1vkRTSKfRe/qCiyExPZduddZOzZS8S993DJ1GcxDIOTH3zAH089hSUwkIvnzKHKNdeQm5XFsjvu5uCXX1Ova2duXfoOPv7+7F/3LYtuv5uczCxumz2DlgP7Y5omf0yYwMl338Xvkkuo9847+IaF4Uw4TNak0ZjHEvG9YzC2u4e5Yv5+NeYXSyGwKpb+ozAiojCzM3F+NA8SD8IlV2C5cSCG5riJSHlTCee4AQwZMuQf57gtW7as2O2W70FkkUI4PlwEaafAkYsj5n+YubnYExNJnPkaOJ3YDx8hcc4bAPzx1kIy9u4D0yRh/gIydu8BIGHSJEy7HUdKCkdffBGAPR/FcfDLrwE48MVX7P3kcwA+f+oZMk+lkJuVxUejnwAge/duTr7zDpgm2fv2uY+xyvkwBvN4IphOcpYuwHkyGdORi7n6PXA6ID0V57qPATB3bXYlbQB7f4Ej+y7Exyci4lmV9OSEoKAgQkJC/vYqCfW4ScVTNRgwXF3tAYFgsWAJCMCw+WLm5IJp4lOjOgC+1au7zxnFYsEn2DX3zVqtmvtQd2tN1xy3gJrV8z0mINR1ykKV0JquYVbTJLCmq8wSFHR29ajTiU91171GUEje8wzw9cHw83fNX/Pzh+xMV53AvPl3/lXyvy//QEREyh2jBIsNjPLfrzRz5kzNcRMpCp9+w8g1TUg9ifWWQRgWCz4hIVzy1psk/m8uAZdeQu2HHgSgzn33YD+ayOlffqXOfffgFxEBQNTcuRx96SUsVatSe9w4V9kN13P9c8+y77OVXNKrJ1HXXwvAra+/zCePP4097TQ9/jMBAFudOkT+97+cePddApo3p8ZA1ypW3zsGY6an4TwSj2+fuzACXcmZ5a5ROL/5GKqGYHRxrUA1Lr0STiVhHt6L0aglRq2zCy1ERKTsKs15yVpVKiIiUkGViVWlL/wfwQHF288sNTObGo+/Vm7zDq0qFRERkfKpEi5OmDdvHtWrV//niiVQ/geRRUREpOwyLCW7yrFhw4Zx8uRJ99f9+vXj6NGjHmm7fH8yIiIiUrZZjJJd5dhfZ6HFxcWRnp7ukbaVuImIiEjpuYA9brNmzaJ+/fr4+/vTqlUr1q1bd866CQkJDBgwgEaNGmGxWBg9enSh9WbMmEGjRo0ICAggMjKSMWPGkJWVVaL4PEGJm4iIiJR7sbGxjB49mnHjxrFlyxY6depEz549iY+PL7R+dnY2tWrVYty4cVx55ZWF1nn33Xd58sknmThxIjt27ODNN98kNjaWsWPH/m0shlHwxBtPrTTV4gQREREpPRdoccLLL79MdHQ0Q4cOBVw9ZStXrmT27NlMmzatQP169erx6quvAjB//vxC29y4cSMdOnRgwIAB7nvuuusuNm3a9LexmKaZ7+SErKwshg8fTpUq+ffnLMnJCUrcREREpPRYSrABb1791NTUfMV+fn6FHiNlt9vZvHkzTz75ZL7y7t27s2HDhuI9+086duzIO++8w6ZNm7jmmmvYv38/cXFx3HPPPX9731+/PzBvL09PUOImIiIipec8etwiIyPzFU+cOJFJkyYVqJ6UlITD4SA8PDxfeXh4OImJicV79p/079+f48eP07FjR0zTJDc3lwcffLBAgvhXb731Vomf+U+UuImIiEjpKclig7z6hw4dyrcB7z8d2v7XeWSmaZ7X3LKvv/6aKVOmMGvWLNq0acPevXsZNWoUERERjB8/vsTtng8lbiIiIlJ6DErQ4+b6T3BwcJFOTggNDcVqtRboXTt27FiBXrjiGD9+PIMGDXLPm7v88stJT09n2LBhjBs3Dktxh4A9QKtKRUREpFyz2Wy0atWK1atX5ytfvXo17du3L3G7GRkZBZIzq9WKaZoF9mq7UNTjJiIiIqXnPBYnFMcjjzzCoEGDaN26Ne3atWPu3LnEx8czfPhwAMaOHcuRI0dYuHCh+56tW7cCcPr0aY4fP87WrVux2Ww0bdoUgN69e/Pyyy/TokUL91Dp+PHjufnmm7FarcWO0ROUuImIiEjpuUDbgfTr14/k5GQmT55MQkICzZs3Jy4ujqioKMC14e5f93Rr0aKF+/XmzZtZvHgxUVFRHDhwAICnn34awzB4+umnOXLkCLVq1aJ3795MmTKl2PF5imEWoa8vNTWVkJAQUlJSijTWLCIiIt7nzX+/zzz7xJzxBAf4F+/ezCxqPPCs8o5CqMdNRERESo9RgrNHPXTKQEWkxE1ERERKz3lsByIF6ZMRERERKSfU4yYiIiKl5wItTqgslLiJiIhI6dFQqUcpcRMREZHSYynB4oTi1q9ElLiJiIhI6dFQqUcpcRMREZHSo6FSj9InIyIiIlJOqMdNRERESo/muHmUEjcREREpPYZRgqFSJW7nosRNRERESo8WJ3iUEjcREREpPVqc4FFK3ERERKT0aI6bRymlFRERESkn1OMmIiIipUdDpR6lxE1ERERKjxYneJQSNxERESk9FovrKu49UiglbiIiIlKKStDjhnrczkWJm4iIiJQezXHzKH0yIiIiIuWEetxERESk9GhxgkcpcRMREZHSo8UJHqXETUREREqPetw8SombiIiIlB7DKMHiBCVu56LETUREREqPetw8SoPIIiIiIuWEetxERESk9GgfN49S4iYiIiKlx2K4ruLeI4VS4iYiIiKlRz1uHqXETUREREqPFid4lBI3ERERKT3qcfMofTIiIiJSIcyaNYv69evj7+9Pq1atWLdu3TnrJiQkMGDAABo1aoTFYmH06NGF1jt16hQjRowgIiICf39/mjRpQlxcXCm9g3+mxE1ERERKjWEYJbqKKzY2ltGjRzNu3Di2bNlCp06d6NmzJ/Hx8YXWz87OplatWowbN44rr7yy0Dp2u51u3bpx4MAB3n//fXbt2sW8efOoW7dusePzFA2VioiISOm5QEOlL7/8MtHR0QwdOhSAGTNmsHLlSmbPns20adMK1K9Xrx6vvvoqAPPnzy+0zfnz53PixAk2bNiAr68vAFFRUcWOzZPU4yYiIiKl50ziVtyrGOx2O5s3b6Z79+75yrt3786GDRtKHPpHH31Eu3btGDFiBOHh4TRv3pypU6ficDhK3Ob5Uo+biIiIlB6jBPu45Q2Vpqam5iv28/PDz8+vQPWkpCQcDgfh4eH5ysPDw0lMTCzes/9k//79rFmzhrvvvpu4uDj27NnDiBEjyM3NZcKECSVu93yox01ERERKz3n0uEVGRhISEuK+ChvyzPeov8yNM02zRPPlznA6nYSFhTF37lxatWpF//79GTduHLNnzy5xm+dLPW4iIiJSJh06dIjg4GD314X1tgGEhoZitVoL9K4dO3asQC9ccURERODr64vVanWXNWnShMTEROx2OzabrcRtl5R63ERERKT0nNmAt7gXEBwcnO86V+Jms9lo1aoVq1evzle+evVq2rdvX+LQO3TowN69e3E6ne6y3bt3ExER4ZWkDZS4iYiISGkyjBIMlRZ/ePORRx7hjTfeYP78+ezYsYMxY8YQHx/P8OHDARg7diyDBw/Od8/WrVvZunUrp0+f5vjx42zdupXt27e7v//ggw+SnJzMqFGj2L17N59++ilTp05lxIgR5/eZnAcNlYqIiEjpuUBHXvXr14/k5GQmT55MQkICzZs3Jy4uzr19R0JCQoE93Vq0aOF+vXnzZhYvXkxUVBQHDhwAXHPsVq1axZgxY7jiiiuoW7cuo0aN4oknnih2fJ5imKZp/lOl1NRUQkJCSElJyTfWLCIiImWXN//9PvPsk9+sILhqleLdezqd6tf2Ud5RCPW4iYiISOmxlGA7kOLWr0Q0x01ERESknFCPm4iIiJSeC3TkVWWhxE1ERERKzwVanFBZKHETERGR0qMeN49S4iYiIiKlRz1uHqXETUREREqPetw8Sp+MiIiISDmhHjcREREpPRaL6yruPVIoJW4iIiJSagzDwCjmnLXi1q9MlLiJiIhI6TlzyHxx75FCKXETERGR0qNVpR6lxE1ERERKUQlWlWrt5DnpkxEREREpJ9TjJiIiIqVHQ6UepcRNRERESo+2A/EoJW4iIiJSetTj5lFK3ERERKT06Mgrj9InIyIiIlJOqMdNRERESo+GSj1KiZuIiIiUIiPvKu49UhglbiIiIlJ61OPmUUrcREREpPQocfMoJW4iIiJSijRU6klaVSoiIiJSTqjHTUREREqPhko9SombiIiIlB6NlHqUEjcREREpRcrcPEmJm4iIiJQeDZV6lBI3ERERKT0GJUjcSiWSCkGrSkVERETKCfW4iYiISCnSHDdPUuImIiIipUdz3DxKQ6UiIiJSiowSXsU3a9Ys6tevj7+/P61atWLdunXnrJuQkMCAAQNo1KgRFouF0aNH/23bS5YswTAM+vTpU6LYPEWJm4iIiJSeMz1uxb2KKTY2ltGjRzNu3Di2bNlCp06d6NmzJ/Hx8YXWz87OplatWowbN44rr7zyb9s+ePAgjz32GJ06dSp2XJ6mxE1ERERKzwVK3F5++WWio6MZOnQoTZo0YcaMGURGRjJ79uxC69erV49XX32VwYMHExIScs52HQ4Hd999N8888wwNGjQodlyepsRNREREyqTU1NR8V3Z2dqH17HY7mzdvpnv37vnKu3fvzoYNG84rhsmTJ1OrVi2io6PPqx1PUeImIiIipajkc9wiIyMJCQlxX9OmTSv0CUlJSTgcDsLDw/OVh4eHk5iYWOLIv/32W958803mzZtX4jY8TatKRUREpNQYhoFRzKHPM/UPHTpEcHCwu9zPz69I951hmmaxn31GWloaAwcOZN68eYSGhpaojdKgxE1ERERKz3lsBxIcHJwvcTuX0NBQrFZrgd61Y8eOFeiFK6p9+/Zx4MABevfu7S5zOp0A+Pj4sGvXLho2bFiits+HhkpFRESkFJX+diA2m41WrVqxevXqfOWrV6+mffv2JYq6cePGbNu2ja1bt7qvm2++mc6dO7N161YiIyNL1O75Uo+biIiIlKKSrBIt/vDmI488wqBBg2jdujXt2rVj7ty5xMfHM3z4cADGjh3LkSNHWLhwofuerVu3AnD69GmOHz/O1q1bsdlsNG3aFH9/f5o3b57vGdWqVQMoUH4hKXETERGRcq9fv34kJyczefJkEhISaN68OXFxcURFRQGuDXf/uqdbixYt3K83b97M4sWLiYqK4sCBAxcy9GIxTNM0/6lSamoqISEhpKSkFGmsWURERLzPm/9+u599YBfBwUHFvDeNkHqNlHcUQj1uIiIiUop0yLwnKXETERGR0qND5j1KiZuIiIiUHnW4eZQSNxGRSiIuLo41a9ZwySWXMGDAAM0dkgtEmZsnaR83EZFK4L333qNXr168/fYCHn74YS6++GI2bdrk7bBEpJiK1eOWmppaWnGIiEgp2rx5MwBBQVUxDDh+PJnNmzfTuHFjL0cmpalM/LutOW4eVaTtQLKysqhfv/55HdQqIiIiF17t2rX5/fff8ff3v6DPdW8Hcnh/ybYDuaiBtgMpRJF63Pz9/fn999+x2+2lHY+IiIh4kM1mu+BJW36a4+ZJRR4q9ff39/L/8SIiIlLuGJRgqLRUIqkQtKpURERESo/muHmUVpWKiIiIlBPqcRMREZFSpDlunqTETUREREpN6unTxR76TD19upSiKf+UuImIiIjH2Ww2ateuTeRlzUp0f+3atbHZbB6Oqvwr0j5uIiIiIsWVlZVV4q3EvL+NSdmkxE1ERESknNCqUhEREZFyQombiIiISDmhxE1ERESknFDiJiIiIlJOKHETERERKSeUuImIiIiUE0rcRERERMqJ/wc8PuoYD4RXkQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 800x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# ftle_array = FTLE.compute_ftle(x0,y0,x1,y1, Td)\n",
    "fig = plt.figure(figsize=(8,4))\n",
    "ax = plt.axes(projection=ccrs.PlateCarree())\n",
    "ax.add_feature(features.LAND)\n",
    "ax.add_feature(features.COASTLINE)\n",
    "plt.scatter(ds['lon'][:,0], ds['lat'][:,0], c=ftle_array, cmap='Reds', s=3)\n",
    "ax.set_xlim(10, 47)\n",
    "ax.set_ylim(-47, -25)\n",
    "cbar = plt.colorbar()\n",
    "cbar.set_label(\"FTLE (1/days)\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Land mass in the release grid**\n",
    "\n",
    "If the region you are interested in also has a land mass, you can do the following:\n",
    "\n",
    "1. Create a regular array of particle locations over the whole area (land and water).\n",
    "2. Use the land mask from the ocean model data to create a boolean matrix corresponding to the released particles and extract the ocean only particles.\n",
    "3. After the simulation, use the boolean mask array to map the released particles initial (and final) locations to an array of the boolean matrix shape. Non-ocean particles are set to NAN.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}