Ft hist and contour vanilla added

notebooks/Ft_computation.ipynb

+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "2f6c6cc8-bc89-4a14-9dea-7ff0bf91d46d",
+   "metadata": {},
+   "source": [
+    "## AGENDA\n",
+    "\n",
+    "In this notebook, we compute the time-domain amplification factors $\\widetilde{F}(t)$.\n",
+    "\n",
+    "\n",
+    "`Author: U. Deka`\n",
+    "\n",
+    "## Table of contents:\n",
+    "- [Histogram method](#Ft_histogram)\n",
+    "- [Contour method](#Ft_contour)\n",
+    "- [Testing $\\widetilde{F}(t)$](#Ft_test)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "ebd281a9-b148-417f-9ac3-bf1279193b6b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import sys\n",
+    "sys.path.append(\"../utils/\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "d1cb591e-bb89-48fe-a44c-a5da34d52895",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "from plotsettings import *\n",
+    "import warnings\n",
+    "import time\n",
+    "from tqdm import tqdm\n",
+    "from cycler import cycler\n",
+    "from matplotlib import colors\n",
+    "from lenses import ChargedLens\n",
+    "import warnings"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "6e9164ad-3ebf-4233-b63d-b2d5429a3b5f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "style = 'Notebook'\n",
+    "plotsettings(style)\n",
+    "plt.rcParams['axes.prop_cycle'] = cycler(color=dark2)\n",
+    "warnings.filterwarnings('ignore')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "fa8766da-6ca6-4ac9-b709-dd9b09026d43",
+   "metadata": {},
+   "source": [
+    "## Histogram method <a class=\"anchor\" id=\"Ft_histogram\"></a>"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "77e7c1a1-3f1a-46b0-bc80-ac01cfee4e9f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def Ftilde_hist(lens, t_max, dt, dx_fac=1, verbose=False, dx1=None, dx2=None):\n",
+    "    \"\"\"\n",
+    "    Computes F(t) with the histogram method.\n",
+    "    Params:\n",
+    "    ---------------\n",
+    "    lens: object of a charged lens class\n",
+    "    t_max: maximum time delay to be resolved (relative to global minima)\n",
+    "    dt: bin width of histogram\n",
+    "    \"\"\"\n",
+    "    #!FIXME we assume that atleast one image exists\n",
+    "    minima_td = min(lens.img_td)\n",
+    "    t_max += minima_td\n",
+    "    x_max = 1.2 * np.sqrt(2 * t_max)  \n",
+    "    dtdx1, dtdx2 = lens.gradTimeDelay(x_max, x_max)\n",
+    "    if not dx1:\n",
+    "        dx1 = dt / dtdx1 * dx_fac\n",
+    "    if not dx2:\n",
+    "        dx2 = dt / dtdx2 * dx_fac\n",
+    "    x1_ = np.arange(-x_max + lens.y1, x_max + lens.y1, dx1)\n",
+    "    x2_ = np.arange(-x_max + lens.y2, x_max + lens.y2, dx2)\n",
+    "    X1, X2 = np.meshgrid(x1_, x2_)\n",
+    "    td_samples = lens.timeDelay(X1, X2).ravel()\n",
+    "    td_samples = td_samples[td_samples >= minima_td]\n",
+    "    td_samples -= np.min(td_samples)\n",
+    "    bins = int(max(td_samples) / dt)\n",
+    "    td_hist, bin_edges = np.histogram(td_samples, bins=bins)\n",
+    "    t_mid = np.zeros_like(bin_edges)\n",
+    "    t_mid[1:] = (bin_edges[1:] + bin_edges[:-1]) / 2.0\n",
+    "    Ft = np.zeros_like(bin_edges)\n",
+    "    Ft[1:] = dx1 * dx2 * td_hist / dt / 2.0 / np.pi\n",
+    "    Ft[0] = np.sqrt(lens.img_mag[np.argmin(lens.img_td)])\n",
+    "    if verbose:\n",
+    "        print(f\"Tmax = {t_max} | dt = {dt} | xmax = {x_max} | dx1 = {dx1} | dx2 = {dx2}\")\n",
+    "        print(f\"Number of td samples = {len(td_samples)}\")\n",
+    "    cut_idxs = t_mid<=t_max\n",
+    "    return t_mid[cut_idxs], Ft[cut_idxs]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "4525fec6-d99d-4802-abb3-a152d0af930e",
+   "metadata": {},
+   "source": [
+    "## Contour method <a class=\"anchor\" id=\"Ft_contour\"></a>"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "6ab60a78-7c19-40d6-ad78-612eb748ee64",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def Ftilde_contour(lens, t_max, dt, dx_fac=1, verbose=False, dx1=None, dx2=None, t_min=0):\n",
+    "    \"\"\" \n",
+    "    Computing the raw F(t) \n",
+    "    for charged lens using contour summation method \n",
+    "    Params:\n",
+    "    ---------------\n",
+    "    lens: object of a charged lens class\n",
+    "    t_max: maximum time delay to be resolved (relative to global minima)\n",
+    "    dt: bin width of histogram\n",
+    "    \"\"\"\n",
+    "    #!FIXME: We assume that atleast one image exists\n",
+    "    minima_td = min(lens.img_td)\n",
+    "    t_max += minima_td\n",
+    "    x_max = 1.2 * np.sqrt(2 * t_max)\n",
+    "    dtdx1, dtdx2 = lens.gradTimeDelay(x_max, x_max)\n",
+    "    if not dx1:\n",
+    "        dx1 = dt / dtdx1 * dx_fac\n",
+    "    if not dx2:\n",
+    "        dx2 = dt / dtdx2 * dx_fac\n",
+    "    x1_ = np.arange(-x_max + lens.y1, x_max + lens.y1, dx1)\n",
+    "    x2_ = np.arange(-x_max + lens.y2, x_max + lens.y2, dx2)\n",
+    "    X1, X2 = np.meshgrid(x1_, x2_)\n",
+    "    n = int((t_max-minima_td) / dt)\n",
+    "    Tarray = np.logspace(0, np.log10(t_max+1), n)\n",
+    "    Tarray -=1\n",
+    "    Tarray = Tarray + t_min + minima_td\n",
+    "    plt.figure()\n",
+    "    contour_levels = plt.contour(X1, X2, lens.timeDelay(X1, X2), Tarray)\n",
+    "    level_coords = contour_levels.allsegs\n",
+    "    Ft = []\n",
+    "    for level in level_coords:\n",
+    "        res_sum = 0\n",
+    "        for i in range(len(level)):\n",
+    "            x_coord, y_coord = np.vstack(level[i]).T\n",
+    "            ds = np.sqrt(np.diff(x_coord)**2 + np.diff(y_coord)**2)\n",
+    "            gtd = lens.gradTimeDelay(x_coord, y_coord)\n",
+    "            gradnorm = np.linalg.norm([gtd[0],gtd[1]], axis = 0)[:-1]\n",
+    "            res = np.divide(ds, gradnorm)\n",
+    "            res_sum += np.sum(res) / 2. / np.pi\n",
+    "        Ft.append(res_sum)\n",
+    "    plt.close()\n",
+    "    Tarray-=minima_td\n",
+    "    Ft[0] = np.sqrt(lens.img_mag[np.argmin(lens.img_td)])\n",
+    "    if verbose:\n",
+    "        print(f\"x_max = {x_max} | t_max = {t_max} | dx1 = {dx1}, dx2 = {dx2}\")\n",
+    "        print(f\"Number of td samples = {len(Tarray)}\")\n",
+    "    return Tarray, np.array(Ft)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ca738bde-3fe1-4e8f-a96d-b69b76ca3086",
+   "metadata": {},
+   "source": [
+    "## Testing $\\widetilde{F}(t)$ <a class=\"anchor\" id=\"Ft_test\"></a>"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "id": "9a43cfbe-15a6-482b-a611-b3af5571132a",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Tmax = 29.97804294466326 | dt = 0.020943951023931952 | xmax = 9.291757835879613 | dx1 = 0.002396783848278275 | dx2 = 0.0022670646756315508\n",
+      "Number of td samples = 63566886\n",
+      "x_max = 9.291757835879613 | t_max = 29.97804294466326 | dx1 = 0.025, dx2 = 0.025\n",
+      "Number of td samples = 1432\n",
+      "Time taken to compute Ft_hist 8.691 s. and time taken to compute Ft_cont 9.664 s\n"
+     ]
+    }
+   ],
+   "source": [
+    "qtest = 0.1\n",
+    "ytest = 0.5\n",
+    "\n",
+    "lens =ChargedLens(qeff=qtest, y1=ytest)\n",
+    "T_MAX =30\n",
+    "DT_FAC = 7.5\n",
+    "W_MAX = 20\n",
+    "dt = np.pi / W_MAX / DT_FAC\n",
+    "\n",
+    "t1 = time.perf_counter()\n",
+    "t_hist, Ft_hist = Ftilde_hist(lens, T_MAX, dt, dx_fac=1, verbose=True)\n",
+    "t2 = time.perf_counter()\n",
+    "t_cont, Ft_cont = Ftilde_contour(lens, T_MAX, dt, dx1=0.025, dx2=0.025, verbose=True)\n",
+    "t3 = time.perf_counter()\n",
+    "\n",
+    "print(\"Time taken to compute Ft_hist %.3f s. and time taken to compute Ft_cont %.3f s\"%(t2-t1, t3-t2))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "id": "21f428d3-2346-4d5c-b94c-0fd503e51900",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk8AAAG6CAYAAAARcKJZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABc4klEQVR4nO3deXhbV50//rckS/Lua2dvVl/Fbdq0aSMn0FKglMpNWRqWSgl8YQYYGqsM68DUamBm2HFktvl9v4VBNjDsTGINS4CWVgqFshUSqWnTJm0dyVmaPZavd8ta7u8PRYpky7ZkyZau/X49j55zdXXuR0e5sfXxOeeeq5JlWQYRERERZURd6AYQERERKQmTJyIiIqIsMHkiIiIiygKTJyIiIqIsMHkiIiIiygKTJyIiIqIsMHkiIiIiygKTJyIiIqIslBS6AfPVxo0bYTAYCt0MIqKZGx0FTp4E1q4FSkvzFwvIX1yiPPP5fHj++eenrMPkKUt+vx+iKE5bz2AwYP/+/XPQIiKiWeL1Ao2NwE9/ChiN+YsF5C8uUZ5t37592jqKHLZzu91obGyESqWCwWBAW1tbRsdZLBaoVKq0D6vVOqG+wWCYUI+IaMFYuRL46ldjZT5j5TMuUQEorufJ6XSitbUVO3fuxJYtW7Bv3z7YbDb09PTAbrdPepwkSXC73Whubp4wnGaz2WCxWCa8jyiKMJvNiX0GgyGjXicionlh2TLg4x+fnVj5iktUACql3Ri4sbERHo8n8VySJNTX10OSJEz1Udrb22EymSYkP16vF42NjROObWpqQmdnJwRBmFE7t2/fzmE7IlK23l7A7QZMJqC2Nn+xgPzFJcqzTL6/FTVs5/V6J/QuCYKAHTt2AIjNR5pMc3Nz2l6jvXv3pvQuxd/H7XbDZrPB6XRCkqTcG09EpDTd3cCOHbEyn7HyGZeoABSVPBmNRpjif7UkifcOzWRIzel0YufOnSn79u7dCyDWW2WxWFBbW4v29vas4kajUfT396c8gsFg1u0jIiKi4qKo5GkyXq93Qu9RJvx+P/x+/4Rj7XY7ZFmGz+dDS0sLAMBqtWaVQHV1daGmpibl0dramnUbiYiIqLgoPnmKJ0AdHR1ZH+twONL2ZMWJogi73Q6PxwNBENJekTeZhoYG9PX1pTx2796ddRuJiIiouCg+ebJarTOe2O10OjNKiIxGY2Ku1VTzqpKp1WpUV1enPPR6fdZtJCIqmLIyYPPmWJnPWPmMS1QAiluqIFl8iQHjDBZZi/dYTdXzlCxej5PHiWjBuP762OKWsxErX3GJCkCxPU9OpxMGgwHNzc0zOj4+ZJdtj9VMEjUiIiKaPxSZPLndbgQCgQmJk9vtTnk+1RCb0+mcsDDmVJxO54wmpRMRKdbTTwN6fazMZ6x8xiUqAMUlT263G1arFZIkoa2tLfGwWq1wOByJem1tbTAYDHA6nRNixIfs4utDJfN6vYlbvsSH6LxeL1wu14wmpRMRKZYsA2NjsTKfsfIZl6gAFDXnyev1oqmpCUBsvtN4nZ2diW1RFCEIQtq1nxwOB4xGY9ohO1EUIYoiWltb0dramridi8vlyt8HISIiIsVSVPJkNBqnvAVLMrPZPOkw21T3wBMEgYkSERERTUpxw3ZERDS3zj68A9HRwUI3g6hoMHkiIqL0rr8ePfZ/wNCYH8HTz+YcC889FyuTt4kUSFHDdkRENIfKyhBZWQf5qApQa3KOhY0brz5P3iZSGPY8ERFReidPoqrDhZJBGSpVjl8XJ08C998fK5O3KS2n04mmpiaoVKpp63q9XtTW1sLLhUfnDJMnIiJKr6cHZX84Ck0Qufc89fQA3/lOrEzeprSyWVcwEAjM6D0yvd0YTcRhOyIiml6uPU+UtXRL7aRjMpnQ29ubdfympib4fL6sjyP2PBERUQZUufY8UVFpampiz1MO2PNERETTK7LkaSQ8huPSpUI3Y0rrhSUoK9HlHEeSJNhsNuzbtw8A0NHRkTKsJ0kS9u3bB1EUEzex9/v9sNls2Lp1Kw4ePAi3253onXI6nYnEyWq1QhCExPqHbrcbnZ2dEAQBXq8XgiCgo6NjwqLS8XoAsG/fPkiSlFiY+sCBAwCA9vZ2OBwOeDwe7Nq1C263GwcOHIAoirDZbBAEAZIkwe/3w263J+4d6/f74XQ64XA40NnZmfgsgUAAu3fvRktLC5xOJ1pbW+H1etHc3Jxyh5G5wOSJiIjSW7YMg/fciEjp87lPGF+2DHjooVgJpG7PwHHpEt7wq/+XW5tm2aP3fhg3LV6ZcxybzQaLxQKbzYampibs2rUrkTx5vV7YbDa43W7Y7fZE8mSxWLB79+5EvfjdOYDYfKq9e/fC7/enJB3t7e2w2+0pQ3lNTU2or69Hd3d3IoGK3+3D5/NBFEVYrVY0NjbCZDIlEiq32w2Hw5GSxPn9fgQCgUQi5PF4AMQSuLvuuiuR3NXV1cHn88Hv96O1tRVWqxUejwc2mw02mw0+nw9NTU04cOAA9u3bB6vViqampjm9/yyTJyIiSm/lSgy8dTPCfz0KZHDV13Sx0Np69Xny9gysF5bg0Xs/nFubZtl6YUle4tjt9kTiYjab0dbWBq/XC6PRCKPRCLvdjsbGxpRjvF4vDh48mEgoprqzBnC1d2v37t0p+x0OBwwGA3bt2pVIjOIJV3xOltFohCiKKcOAJpMp0VaLxQKTyYSWlpbE8fFeJgAwGAyQJCnReyUIApqamtDe3g6r1ZqSELa3t0MQhMTn2rFjB6xWK1wuF5MnIiIqAgMD0L14FqqQnPtNfAcGAI8HiH/Jx7erqmYUrqxEl5deHSVIHjLbunUrgNQr7NLdpzWeuIwfEpvMoUOHIEnShEnq8fu9ut3uCe8XT3aAWAIlSVLa2Fu2bEl5nnwfWqfTmfEt0eJtW7Ro0YS2zPSKw5nihHEiIkqvqwuLvn4Aun5ARo7JU1cXcOedsTJ5m2ZFZ2dnYm6QwWBAe3v7lPXjvUbpkhBRFFMSo/g8KZvNltgXHzbMVLxXSRRFWCyWjI8rFux5IiKi6eXa80Rzzm63o6mpCVarFVarFXV1dZMObcV7dSZbuiC5d0sURXg8nsQ8LADweDwZL63Q1NQEQRASPVCHDh3K9CMVDfY8ERFRBpg8KUlbWxuA2Nyj+MTuvXv3TlrfZDJBEIS0PVR+vx/Nzc2J5/H5UR6PB3a7HXa7PePEyev1wu12J4Yf4/GAuR96ywWTJyIimh57nuZcPJlInogd304eRotv9ySt2D7+0n1BEFISlniy43Q64fV64ff70dHRAUmSUhIop9MJIHXCeXt7O9xuN2w2G9ra2hKP5HlR49uY3I54+7xeL9rb2xNzntxud+IWM1MlUul6xyabbzVrZJoV9957b6GbQESUm2eflcM1pXL3dpU8eurZnGPJK1fGyuRtSquzs1MWBEEGIJvNZtnn88kul0sWRVEGIBuNxsQ+k8kkA5BFUZQdDocsy3Liud1ul5ubm+WWlpaU+D6fTxZFURYEIXGMLMuyx+ORTSaTbDab5ZaWlgnHybIs9/b2ykajUTYajYk2xh/x+g6HI6X9LpcrcXxLS0viM3R2dsq9vb2yKIopn8loNMoAZJPJJHs8Htnn88lms1kGIAuCkDiuubk5ZV8+ZPL9rZJl/jkxG7Zv3479+/cXuhlERDl5+WtvwvCzv8Xazz8N/epNhW4OFYG2tjYYjcbEEgJxfr8/sWyAkmXy/c0J40RENLloNFby72zC1UU50/W7iKKYshjnfMY5T0RElN6RI1j+hSeg65XTfllmGwurVsXK5G1SlPhcJIvFkrjNS/x2KlarNWVi+XzG5ImIiNILhVDSF4QqCuR8tV0oBJw5EyuTt0lR4rdg8fv9sFgsaGxshNVqRSAQgMPhSLtg53zEYTsiIpoeh+3oCrPZPKe3QilG7HkiIqIMMHkiimPyRERE02PPE1ECh+2IiCi9hgZcsm7B2IAn9+SpoQF44olYCaRuEykMkyciIkqvqgpBQx3k51W53xi4qgp43euuPk/eJlIYDtsREVF6Z86g+tGXUDIk597zdOYMsHt3rEzeJlIgJk9ERJTehQuofuIENKPIPXm6cAHYsydWJm8TKRCTJyIiygAnjBPFMXkiIqLp8Wo7ogQmT0RElAEmT0RxvNqOiIjSW7QIQ1uvQUR/Lveep0WLgPe/P1YCqdtECsPkiYiI0lu7Fr3mjQgfPZ/7jYHXrgW+/e2rz5O3iRSGw3ZERJTeyAg05wegCsvIedhuZAR4/vlYmbxNRcfv9xe6CUWPyRMREaV37BhWfP0p6PqQ+7DdsWPAjTfGyuRtKjpNTU2FbkLRY/J0BTNtIqIp8Gq7BaGpqYnfhxlQ5Jwnt9sNm80Gr9cLURRhtVrR0tKS8fEGg2HCfw6fz5fy3O/3w2azQRRFSJIEQRBgt9vz0n4iIuUpruQpGhzG2LkXCt2MKelWbIBaX55TDK/XC4fDAQA4dOgQTCbThO8it9uNzs5OCIIAr9cLQRDQ0dEBQRAAxL7PHA4HnE4nPB4PWltb0d7ejrq6OjgcDphMJgCA0+lMfDdardaU773p3qO9vR02mw2SJCXmx7W1taG1tRWSJKG3tzdRz+FwwOPxYNeuXXC73Thw4ACMRmNO/05zTXHJk9PpRGtrK3bu3IktW7Zg3759sNls6OnpySi5cTqdEEURZrM5sc9gMEAUxcRzv9+PxsZGdHR0JOpZLBZYLBZ0dnbm/0MRERW7Iut5Gjv3Ak59ZmuhmzGlNZ85iNJ1M08K3G437HY7XC5X4nl8SC3+fdfe3g673Z7SAdDU1IT6+np0d3dDEATU1dXB6/XC7/dj165d2LlzJ3bu3Jn4XosnNmazGXv37k0kW3GZvEdzczNcLhecTmeiTktLC3w+H9rb2wHEkj+Hw5HonNi6dSv8fj8CgcCM/40KRXHJU2trKzweT+K53W5HfX092traMkqeHA5HInuejM1mQ11dXUqCZbfbYTAY4HQ6U/YTEc1bKhWiGtWVJzkmTyoVoNPFSiB1ewZ0KzZgzWcO5tamWaZbsWHGx0qSBIvFggMHDiT2mUwmNDc3Y+vWrYk6NpsNu3fvTjnW4XDAYDBg165die+7pqYmuN1u7N69O9HLYzab0dbWBr/fn9KBML4dmbwHANTV1U04Pvm71mQyJd7TYrHAZDJlNWpUTBSVPHm93gkJkiAI2LFjB9rb26f8DxA/Pj7k19TUBJPJNCGJkiQJTqcTzc3NKftFUYQgCHA4HEyeiGhh2LwZf265CcvOPpf7UgWbNwPB4NXnydszoNaX59SrU+wOHToESZImDGcl9wjF64z/3hNFEaIowu12T4ib/J1nMBgAxL73pmtHNu+RiS1btszouGKhqAnjRqMxMTabLP6fYarECQD27t0LINYFabFYUFtbm+hOjIv/R4j/p0qWy38UIiIlSqRMRTZsN995vd5p68TnJ6Ub9orP183VXLyHEikqeZqM1+vNqDfIbrdDlmX4fL5EV6HVak1JoOL/UdIN68W7JDP5zxKNRtHf35/yCOb4lxYR0Zw6dgybv/sSdFIe1nk6dgwwGq8uVRDfprTi30FT/cEe7zAYf8HT+Bi5mIv3UCLFJ09+vx9+vx8dHR0ZHyOKIux2OzweDwRBgNVqnVBnqrHbTC7j7OrqQk1NTcqjtbU14zYSERXcyAiqLoxCFUHuPU8jI8DTT19dJDO+TWnFR1lsNlvK/vgcpHgdQRAmjKAAse+p8dNPZtqObN8juYNhvvZMKT55slqt004An4zRaEzMoYonRPEsO10XZfw/QSbv1dDQgL6+vpTH+Al3RESKwWG7ORW/Ktzr9cJgMMBms8FqteKuu+5K+S7p6OiAJEkpyU38irfkOcI9PT0AUpOZdL1J8e9Ap9OZuEIv0/dobGwEgMQSBDabLfHdGj92viRTik6ebDYbLBZLTutDxLP7+AmN/8dJd4LjCdV0c6sAQK1Wo7q6OuWh1+tn3E4iokKS5Wihm7DgdHZ2JqaYtLe3IxAITOgsMJvN8Hg86OzshMVigc1mw8GDB1MSI7fbnUh2Wltb4ff7U/YlJzlWqxWiKGLXrl04dOhQIomb7j0AoLm5GWazGU6nE1arFU1NTbBYLDCbzRBFEe3t7di3bx+AqwmWUqnknC+hKAyn04lAIJBzt6Tf74fBYEhcSSJJEmpra9Hc3JxyVQMA1NbWQhTFlKUSJrN9+3bs378/p7YRERWU1ws0NuLkm1Wo/fx+VN/y5pxjIf77M76tsMURaf7L5PtbkT1Pbrc7beI0PovNZG7S+HWbBEGA2WxOG0uSJA69EdHCUV+PI9tXIVQJRKI59jzV1wP79sXK5G0iBVJc8uR2u2G1WiFJEtra2hIPq9Wa0lPU1taWWNQSQGLcuK2tLTEk5/V64XK5Jkw2t9vtCAQCKQmUzWZLLPBFRLQg1Nbi0oYqRPUqRKKRnGPBYomVydtECqS4RTLjS9OPvwIBQMqtU+KLWsbnJ8UX9GptbUVrayuam5thMBgSy94niw/N2Ww2uFyuxOKbvLcdES0oFy5g1cEeaJbIiOY6w+PCBeDHPwbe9a7Y8/j2smW5t5Nojil2zlOx45wnIlK8pDlP+n/bi+WvtOQci3OeqNjN2zlPREQ0tyJyjsN2RPMIkyciIppWNMpBCqI4Jk9ERDStKHueiBIUNWGciIjmUE0NLooViGqHIed6tV1NDXDvvbESSN0mUhgmT0RElJ7BgMNvuwbrLvugznXYzmAAkifh8oIaUjAO2xERUXqhEEqGw0BURjTX27OEQsClS7EyeZtIgZg8ERFRekeO4PX/1Q19L3JfJPPIEWDp0liZvE2kQEyeiIhoWrwxMNFVTJ6IiGhaOQ/bEc0jTJ6IiGhaOd8YmGgeYfJERETT4jpPRFdxqQIiIkrv5pvx+AfWYt3gKURznTB+881AXx9QURF7nrxNpDBMnoiIKD2NBiG9BhhW5b5IpkYDVFdffZ68TaQwHLYjIqL0urqw9efnoO2Xc1+qoKsL2LYtViZvEykQkyciIkpvYABLT45AHQLkSDjnWHj88ViZvE2kQEyeiIhoWjnPeSKaR5g8ERHRtKK59jwRzSNMnoiIaFo5Txgnmkd4tR0REaW3ejWO3CGgvKIPcjTHnqfVq4GHH46VQOo2kcIweSIiovSWLMHpmyqwfqA/9zlPS5YAH/zg1efJ20QKw2E7IiJKLxDANS8OQh2Ucx+2CwSAH/0oViZvEykQkyciIkrvxAnc4uqDdjAPc55OnAD+4R9iZfI2kQIxeSIiomnlPOeJaB5h8kRERNOKRqOFbgJR0WDyREREaUXlqwkTe56IruLVdkRElFakvBwDyzSIlkRzn/NUUQHcemusBFK3iRSGyRMREaUVaWjA02+txjWjfUCuydN11wF//evV58nbRArDYTsiIkorKstQyzIAQOacJ6IEJk9ERJRW1OvBa9sl6Hvk3Oc8eb2AShUrk7eJFIjJExERpZUyYVxmzxNRHJMnIiJKK3plyC72hDcGJopj8kRERGlFkuY5yREmT0RxTJ6IiCitKJJ6nmQmT0RxXKqAiIjSCl53LZ6z6KHXj+V+td0NNwBdXcCqVbHnydtECsPkiYiI0orodQhVqaCLqHKf81RaCqxff/V58jaRwnDYLkt+v7/QTSAimhMqfzfW/X4MJQNy7sN23d3Au98dK5O3iRRIkcmT2+1GY2MjVCoVDAYD2traZuVYg8EAlUqV8iAiWjAkCbX+KDRjeViqoLcX+PGPY2XyNpECKW7Yzul0orW1FTt37sSWLVuwb98+2Gw29PT0wG635+1Yp9MJURRhNpsT+wwGA0RRnJXPRURUbCLJCROXKiBKUMly8kIexa+xsREejyfxXJIk1NfXQ5IkTPdRsjm2qakJnZ2dEARhRu3cvn079u/fP6NjiYiKgf+JxyC+/h6cfLMKR7bejjf/xx9nHszrBRobgfjv4Pi20ZifxhLlSSbf34oatvN6vRN6iARBwI4dOwBMPR8pm2O9Xi/cbjdsNhucTickScrTJyAiUo5I5OotWbjCONFVihq2M07yF0q8d2iqIbVsjt27dy8AoL29He3t7QAAh8OB5ubmjNsajUbR39+fsk+v10Ov12ccg4iokEaXLELPzUC4DECuSxWsWAF8+tOxEkjdJlIYRfU8Tcbr9abMTcr1WLvdDlmW4fP50NLSAgCwWq2JRCoTXV1dqKmpSXm0trbOqI1ERIUwtkhAzy1qRMpVUOV6td2KFcBnPhMrk7eJFEjxyZPf74ff70dHR0fejxVFEXa7HR6PB4IgwGq1Zhy7oaEBfX19KY/du3dn3UYiokKRpV6Un5GhHpNz73nq7wceeyxWJm8TKZDikyer1Trjid2ZHms0GhPzpTJd50mtVqO6ujrlwSE7IlISnc+PVW4Z2gHkvs7T8ePAPffEyuRtIgVSdPJks9lgsVgmnc+Uz2NNJhMAcPI4ES0Y0SsTxqNQAcq6MJtoVik2eXI6nTAYDFlN4s7HsTNJ1IiIlCgaDQEAImoNVFzniShBkcmT2+1GIBCYkPy43e6U5+mG2DI9djyn0znjSelEREqU6HlSawAuVUCUoKilCoBYkmO1WmG1WlNureLz+RAIBBLDa21tbbDZbOjs7EwkPZkc6/V6YbFYYLVa0dzcDEEQ4PV64XK50NnZObcfloiogCIaFcaqgJBWB3WuyZNeDxgMsRJI3SZSGEUlT16vF01NTQBic5bGS05uRFGEIAiJ9ZsyPVYURYiiiNbWVrS2tqK5uRkGgwEulyvvn4eIqJgNrluJE29XY7hCn/vVdhs3pk4Q52RxUjBFJU9Go3HaW7DEmc3mlGG2TI8VBIGJEhERgGgkNucpqtFCxWE7ogRFznkiIqLZV/5CF8T/iUIfkHNPnp59FliyJFYmbxMpkKJ6noiIaA6FxlASBKAqyX3CeDgMXL4cK4HUbSKFYc8TERGlJUevXG2nKeGwHVESJk9ERJRWPHmCRgtVrhPGieYRJk9ERJSWfGWdJ1mjzX2pAqJ5hHOeiIgorYEVi3DqDSoML62BWrqcW7BrrwX+8pdYCaRuEykMkyciIkorrFNjdKkKocoKlARyvD1LZSVw221XnydvEykMh+2IiCgt3YXLWHIwCt2QDE2u97Z7+WXg4x+PlcnbRArE5ImIiNLS90qoPQqUjMjQyDkmTxcvAl//eqxM3iZSIA7bERFRepHYJHGVVg81r7YjSmDPExERpZVYqqCkNPeeJ6J5hMkTERGldyV5Umn10ESjiHK5AiIATJ6IiGgSo+VaSNcB4ToBWjmCUC5Dd4sXA//8z7EyeZtIgTjniYiI0hqpq8DZ20oQXr4UJXIUwXAIes0MvzbWrAG+8Y2rz5O3iRSGPU9ERJSWZngYuh6g5Mp0p1A4OPNgw8OA1xsrk7eJFIjJExERpVV5pgfrfxVG5TkJABAcG5l5sBdeABobY2XyNpECMXkiIqL0rkwY15RoAQChUA49T0TzCJMnIiJKS51InvQAgLExDrMRAUyeiIhoMuEQgOSep9FCtoaoaDB5IiKi9KJhhLUqaLQ6ADkO26nVQFVVrEzeJlIgLlVARERpSctK8ecPrMI1N2wAnsgxebrlFqC//+rz5G0ihWHaT0RE6YWDiGq0KLnS8xTmhHEiAEyeiIhoElVnJWz9wcsoO3UBQI5zno4eBTZujJXJ20QKxGE7IiJKqyQYRGXPGIJXFsnMqedpdDSWLI1eScCSt4kUhj1PRESUXjR2tV3JlaUKwrmsME40jzB5IiKitNSR2DpPJSWc80SUjMkTERGlFU+etFeSp0h4rJDNISoanPNERERpjVQBh//xFmxcvx4AEMllwrgoAr/8ZawEUreJFIbJExERpSVro7i8aTXUi5cBAKK5JE+CAGzffvV58jaRwnDYjoiI0iobCGLdE91Q9fQByDF5On8eaG2NlcnbRArE5ImIiNIqHwhh/W+eg+pSDwBAziV5OnsW+OQnY2XyNpECMXkiIqIJItEoSuTYAk8qTQnCKnVuPU9E8wiTJyIimiAYCaMkGkk8D6tLcut5IppHmDwREdEEwUgoJXkKabQA13kiAsCr7YiIKI2RSBgl2ggCt96AOkFAWKMFcllhXBAAszlWAqnbRArD5ImIiCYYCQ5DrpRx4lPvRZ0oIqzRQhXOcZ2nzs6rz5O3iRRGkcmT2+2GzWaD1+uFKIqwWq1oaWnJ6Fi/3w+bzQZRFCFJEgRBgN1un3E9IqL5KDjSj5KIjLKBEDA2hohGC4RyWGF8bAy4eBFYujT2PL6t0+WnwURzSHFznpxOJ2w2G3bu3Inm5mYEAgHYbDbYbLZpj/X7/WhsbMTOnTtht9vhcDjg9/thsVhmVI+IaL4Kjg5ALwHX/59PAc89h4hGC1Ukh+TpueeA1atjZfI2kQIpLnlqbW2Fx+NBS0sLHA4Huru7IQgC2trapj3WZrOhrq4OZrM5sc9ut8PpdMLpdGZdj4hovhobHUx5HtXooM5lzhPRPKKo5Mnr9U4YOhMEATt27AAQ6zGajCRJcDqdMJlMKftFUYQgCHA4HFnVIyKaz8YnT5ESHVS8MTARAIXNeTIajWn3C1eu2BCnuMmk2+0GABgMhgmviaKYeD3TetOJRqPo7+9P2afX66HX6zM6noiokEKjAynPoyU6aNjzRARAYT1Pk/F6vSlDbOnEe6WENJfG1tXVAYj1OmVabzpdXV2oqalJebS2tk57HBFRMQiNDqU8l0v00IRDBWoNUXFRVM9TOn6/H36/Hx6PJ6P68QQoWTxRSh72m67eZL1gcQ0NDTh48GDKPvY6EZFShINDCNYBY2ePQ7esHvKjOmhymTB+yy3A6Cig1caeJ28TKYzikyer1YrOzs60PUXJ4kN6gUBgwmvxniRBEDKuNx21Wo3q6upp6xERFaPI2DCgUkFTXQeo1ZBL9CiJ5NDzpFYDyX9A8o9JUjBFD9vZbDZYLJZpe4GAq8lTuiG3eKIkimLG9YiI5rNocAjaPhnqN74FeOml3JOnl14CXve6WJm8TaRAiu15cjqdMBgMaG5uzqh+POHx+XwTXksehsu0HhHRfBYJDgFhFVRP/hEYHAS0epREc0ieBgeBP/whVgKp20QKo8ieJ7fbjUAgMCFxGn8lXPIcJkEQYDab09aRJAm7d+/Oqh4R0bwWHMSY5uqcJJVWj5JIuIANIioeikue3G43rFYrJElCW1tb4mG1WlPWYGpra4PBYEhZ1NJutyMQCKQkRjabDSaTacKCmJnUIyKat0aHMKZJunVKSSm0uQzbEc0jihq283q9aGpqAoC0t2PpTLrRZHxRy+T5SaIowuPxwGazweVywe/3QxTFCQtvZlqPiGi+Uo0NIVRyNXlS68qhjbLniQhQWPJkNBohy3JGdc1mc9peIlEUU5KsyWRaj4hoPlIHhzBYWwF07AHWrEHJi2XQR0IIRyIo0WiyD7hmDdDRESuB1G0ihZmT5Km/v5+X7RMRKYh6bBjDNRXA/fcDALSlVdBAxvDoIKorarIPuHhxIhaA1G0ihcn7nKf+/n585StfwbZt27Bo0SJoNBrU1tZCo9FAo9Fg27Zt+OQnP4kTJ07k+62JiChP1GNDUId1wLe/DVy+DF1Z7A/goWFpZgEvX07EStkmUqC8JU/9/f34wAc+gNraWrS0tMDlcqG3txeyLKc8XC4X9uzZA4PBgHvuuQcnT57MVxOIiChPNGMjKBvRALt2AadOQVceS56GB3tnFvDUqUSslG0iBcrLsN3PfvYz3H+lC/a+++5DU1MTtmzZAkEQUm5zEggEEvePe/zxx9HZ2QlRFNHe3o73v//9+WgKERHlgXZsBCrtssRz/ZWep5ERqUAtIioeOSdPX/7yl+FwONDR0YH77rtvyro1NbFx8s2bN+O+++6Dw+GA2+3GAw88AJ/Phy996Uu5NoeIiPJAFxqFuqo88by0rAYygNGhvsI1iqhI5DRs99BDDyEQCOD48ePTJk6TMZlMOH78OKLRKL7yla/k0hwiIsoDWZZRGh6FurQysa+sohYAEBzpL1SziIrGjJOnAwcOAABaW1vz0pA9e/bg8uXLOHz4cF7iERHRzIyEQ6gIjwK1i4E77gAqK1F25Qq7sZkmT5WViVgp20QKNONhuy1btuCuu+7KZ1uwZ88e9PWxS5iIqJD6hwLQylFEG9YDv4/duaFisBc9AEKjM0yerr0W+P3vrz5P3iZSmBn3PMXnL6WTS+/RVHGJiGj29fddBACUltcBwSAQjUJ/5Wq78MgMb+YbjSZipWwTKdCs3Ntu165dsxGWiIjmwFB/LHmqPTcAlJYChw9DpdYgqC5BeHRgZkEPH07EStkmUqC8rTD+wAMPwOPxwGQyQZKkKesePnwYt9xyS77emoiI8mho4DIqcHWSeFywRIdocKgwjSIqInnrefrWt76FdevWoa2tDX6/HxqNBvfccw+++tWvTlhNvKenBz/72c/y9dZERJRHowOxlb8rK1OTp7ESPeTgDIftiOaRvA7bdXZ24pZbboEsy6ipqcHjjz+OBx98EAaDAYsWLcI73vEOfOc734HBYEAgEODSBERERSg4GAAAlFXUpewPleghjw0XoklERSXvNwbu7OzEjh07cOjQIXR3d8PtduPxxx/HgQMHsG/fPnR2dibqCoKAf/3Xf813E4iIKAfhoQDCKjXUurKU/aGSUoDDdkT5T55EUYTdbgcA1NfXY9euXYkJ5E8//XRKMrVly5Z8vz0REeUoPNSHEW0ZVDfdBJw+DSxdCgCI6kqhHhuZWdAbb0yJlbJNpDB5T54ATLr+0+bNm7F582Y8+OCDs/G2RESUB/KIhFFtGaDTAatWJfaHdZUoGZthz9O4WCnbRAoz4zlPs7USOFcYJyIqLNVIP8Z05YDfD1gssRJARF8B7UyH7ZJjjYtLpDQzTp7q6+uxc+fOfLYF27ZtgyiKeY1JRETZKRnpw1hZDSBJgNMZKwHIpVXQh2Y4bJcca1xcIqXJaYXxXbt2YevWrTh58mROjThx4gS2bduGlpYWVFdX5xSLiIhyox/pQ7hcmPhCaTVKZ5o8Ec0jOc15MplMAGJzmR544AE0Nzdj3bp1GR9/4sQJ7NmzB52dnThw4AAXziQiKgKlo/0Yrtw08YXyGpSFRue+QURFJucJ4yaTCYcOHcKOHTtgt9shCAJMJhNEUcSiRYsgCALq6uoQCAQgSRJ6enrg9Xpx6NAhSJIEk8kEv9/Pe9oRERUBWZZRGRxCsHrilXDacgFlkTFEwyGoS7QFaB1RcZg2ecrkViqiKOLQoUNob2+H3W5PrOWkUqkm1JVlGQBgNBqxe/du3HfffTNoNhERzYah0SFUh0cwWLMMuOYa4EtfipUA9JWxRTMHBi6hpvaa7AKPi5WyTaQw0yZPX/rSl/DOd74Tb3vb26YN1tzcjObm5sR6TgcPHoQkSQgEAqirq4MgCNi6dStMJhM2b96clw9ARET5czlwCgBQLqwAli8Hdu9OvFZWtQgA0Nd3MfvkaVyslG0ihZk2edqzZw+2bt2aUfIUF1/PiYiIlKX38suoAFBVuzJ2NdyTTwKvfS0gCKioXAwAGOi7mH3g5FhASlwipZn2ajtRFNHb24vf/e53c9EeIiIqoIHeMwCA2sVrYuswveUtifWYKquXAACGBi5lHzg51ri4REqT0VIFgiDAP8l/8kWLFmHbtm3YvXs3F7gkIlK4IekcAKBu0eoJr9UIywAAw4M9c9omomKTUfIkSRKkSRYz6+3tBQDs2LGDV8wRESncmHQeIZUG2jTrPNVUx5Kn4FBgjltFVFymnfPU19c35esqlQqdnZ1c3JKIaB4ID1zCoL4y7dXSJbpSBNUlGBvsLUDLiIrHtD1Phw4dAoApb5syPnHi8B0RkTLJgz0YKa2KPSktBW64IVZeMaItQ3hYyj5wcqw0cYmUZNqeJ5fLBZVKBaPRmPb1dEN1Bw8ehMVigcViwY4dO7hyOBGRQugGLmH0ylV1uOEG4PnnU14PassgD089IpHW+Fjj4hIpybTJk9frhSAIk952pa6ubsK+Xbt2oba2Fjt27IDD4YDJZMLOnTvx9re/PecGExHR7Ckf6sHIyo2Tvh7UVwIz6XkimkemHbbz+/2Je9ilk25cHADMZjMAoLu7G3v37mXiRESkANUjElQ1y2NPDh8Gqqtj5RVj5TXQjkjZB06OlSYukZJklDxNNd/J5/PhiSeeSPtabW0tJ5ITESnE6OgQhLEh6Betie2IRoGBgVh5RaRMgH6kP/vgybHSxCVSkmmTJ0EQYDAYpqxjMpmwaNEifOADH8DPfvazvDWOiIjmzrlzLwEAqpdM/gczKutQFhyYoxYRFaeMVhj3+XxT1pFlGb29vWhvb4fFYoFGo8G2bdsAACdPnpz0OF6VR0RUPC6ejyVPi5evn7SOunIxKsaG5qpJREVp2uTJZDJNuro4EOuZstvtaG5uRn19PWRZhizLcLlckCQJoiim9Er191/t7nU4HPn5FHNgqn8DIqL5oP9i7Pfc8ms2TFpHV70UleEggsHhuWoWUdGZNnmyWq3weDyTvr5jxw48+OCD+Na3voXjx4/D5/MlrrCLJ1LJvVK1tbXYunUrdu7cifb29qwbLEkSbDYbmpqaMj7GYrFApVKlfVit1gn1DQbDhHpERPPd8OVTGNHoUFG1KLZjwwbA44mVV5QJscnkPYGXswueHCtNXCIlmXapgvr6ejgcDnz729/G/fffP+F1m802of6uXbuwa9cuAMCBAwfgcrngdDoTvTcejwderzfrxrrdbrhcLrS1tU267tR4kiTB7Xajubl5wtwtm80Gi8WSss/pdEIUxcTVgkAsmZpq0jwR0XwQ7j2D/jLh6o7ycmDc79rK2hUAAKnnNK5ZcW3mwcfHyvB3OFExmjZ5AoC77rpr0tfq6+unPfauu+7Cnj170NfXh3379iWSqWx7dEwmE0wmU1Y9Vvv27YPH45mQ/MSTt/HLMDgcDnR2dkIQhKzaRkSkdFrpDIaql17dceoUYLcDNhuwJnYFnlC7EmMA+nvPZhc8ORYwIS6RkmR0Y+B8qampwa5du7Bv3z4EAoFpE6/JpFuYczLNzc1pe4327t2b0rsExBIqt9sNm80Gp9M56c2QiYjmo4r+CwgLK6/uuHwZ+OY3Y+UVdYvXAgCGpXPZBU+OlSYukZLMafKUTBCECUN+c8npdGLnzp0p+/bu3QsAKfOzZjIvCwCi0Sj6+/tTHsFgMOd2ExHNhmg0ikVDPVAvWjtlvcrKRQip1BjtuzBHLSMqPgVLngAk5kXNNb/fD7/fP6HnyW63Q5Zl+Hw+tLS0AIhNmJ9JAtXV1YWampqUR2tra17aT0SUb1LvGZRHxlC2bPJlCgBArVZjUFeB8AB7jWjhKmjyVCjxqwEnI4oi7HY7PB4PBEFIe0XedBoaGtDX15fy2L17dy7NJiKaNS+ffBYAsGTVDdPWHdZXITrYM9tNIipaCzJ5cjqdGSVERqMRdrsdQPbrPKnValRXV6c89Hr9jNpLRDTbLpw5CgBYs+6WqzuXLgX+5V9iZZJgWTVUQ1kmT8mxJolLpBQZXW03n8SH7KbqeUoWr8fJ40Q0nw2efwnDJXpUxW8KDACrVgFf+9qEuqGKWuiyHbYbHytNXCKlWHA9T/Ehu2yXIsh0XSkiIiWKXupGX+WS1CVkBgeBv/41VibXrVqCiuHe7N4gOdYkcYmUYl4mT1MNsTmdzgkLY07F6XROmFhORDTflAZOYaRu3JpLL70EvOpVsTKJRrgG1aP9yEpyrEniEimFIpOnQCAw6WttbW0wGAxwOp0TXosP2e3YsWPCa16vFwaDAW1tbYkhOq/XC5fLhY6Ojry1nYio2MiyjEX956Ga5kq7uPLaVSiPjGGAV9zRAqWo5MntdsNqtUKSJHi9XrS1tU24zYsoihAEIe3CmA6HA0ajMe2QnSiKEEURra2tqK+vh81mw6FDh+ByubjaOBHNa+cDL6NubAiVKzdmVL96yToAwIXzx2exVUTFS1ETxuO3Z3E4HJPWMZvNkw6zxa+cS0cQBLhcrpzbSESkNKd8B1ELYOm6zRnVX7RURAhAzwU/1jfcOqttIypGiup5IiKi/Lvc7UUUwBrD1tQXSkqAxYtjZZJlKxoAAP2XT2T+JsmxJolLpBT8n0tEtMAFX34OlysWY0NZdeoLmzYBly5NqF9eWYcRjQ4jgdOZv8n4WGniEikFe56IiBa40gsvoX9xdjdqHyitQTjbmwMTzRNMnoiIFrBoNIolvachr9gw8cXnnwfWr4+V44xU1EHVfzHzN0qONUVcIiXgsB0R0QJ28YIfNaFhjKSbLB4MAj5frBwnUrUYur7zmb/R+FiTxCVSAvY8EREtYKde+hMAYGXDq7I7UFiBqqHJ19wjms+YPBERLWDS8b9isESPtWszW6YgTrdoLepG+zEaGpullhEVLyZPREQL2QkPXq6rR0mWywZULDVAJ0dw4QIXyqSFh8kTEdECJcsyFp1/EaE1N6evsH498NvfxspxaldcCwC4dOZYZm+WHGuKuERKwAnjREQL1PmTz6AqNIza616bvkJ1NbBtW9qXlq+6AecB9Gfa8zQ+1iRxiZSAPU9ERAvU8cOPAgA23PKG9BXOnQM+85lYOU51zTIMafQYvdSd2Zslx5oiLpESMHkiIlqgBrr+hDOVS3DNkrXpK5w7B3z2s5MmOb3ltYj0ZLjKeHKsaeISFTsmT0REC1T5y0cgrbhhxscPVS1FSd/ZPLaISBmYPBERLUDDgwEs6zsLvfjKGccICytQls0q40TzBJMnIqIF6MVnHoMGMtbc1DTjGLolIuqGehCNRvPYMqLix+SJiGgBOn/sCQxp9Lhuw2smr1RbC7zrXbEyjeqVN6A0GsL58y9N/4bJsaaJS1TsuFQBEdECJHcfwrnFIjaXaCevVF8P/OhHk768fO3NkAGc6X4a11yT5sbCU8WaIi5RsWPPExHRAhONRLD0YhfCa6a5JcvoKHD8eKxMY826zYhAhd7TR6Z/0+RY08QlKnZMnoiIFpiXnj+A6tAwFt88zUKVR48CDQ2xMg29vhw95bUYzWTYLjnWNHGJih2TJyKiBcb31F4Ma3TYvOXtOcfqr7kGmssZLpRJNE8weSIiWmBKXvwDXl55E8pKy3OOFVq0FhUS13qihYXJExHRAnLx0kms6TmB0hvvzks83fIGLBm6jHA4lJd4RErA5ImIaAF55nftAIBNd7w/L/GEVTdCK0fxciaTxonmCS5VQES0gIS8v8TJpdfi+mX101c2GgFZnrLKirW3YAzAuZOHsa7emHmsaeISFTP2PBERLRDnL55A/YUXoNm8PW8xV62+EWNqDST2PNECwuSJiGiBOHTAAQBoND2Q2QEvvgjcdlusnERJiRYXKpchfOZY5rEyiEtUzJg8EREtEPLh/Xh52XWoW7IuswOGhoCnnoqVU1VbUg/9peOZx8owLlGxYvJERLQAnDjfhfUXXoTe+Na8x9ZccwOWSmcQiUbyHpuoGDF5IiJaALwHOqCCjM13WfMeu3adEeWRMbx8+rm8xyYqRkyeiIgWAM3h/Ti7/HpULV6T99jrrn0VAODkS3/Je2yiYsTkiYhonjv2wp9x/aUulL3q3dkduG4d8MMfxsopLF95A4Y1Okgnns4sVoZxiYoV13kiIprn/L/9KpZqy7H57g9ld2BdHfDu6RMutVqNy8IqRM9OcaPf8bEyiEtUrNjzREQ0j0VGBnDNc4+he+M2lJZWZnfwpUvAN74RK6cxuvxa1FzyZRYri7hExYjJExHRPPa862GUhoNYs+2j2R98+jTwoQ/Fymno1zVixeBFDA4Epo+VRVyiYrTgkye/31/oJhARzQpZljH0hAPeZRuw9brbZ/W9ll/3aqgBHD/6u1l9H6JioLjkSZIk2Gw2NDU1ZX2swWCASqVKeYzn9/thsVhgs9lgtVphs9ny0Wwiojn3zFOdWNR7GuWvs0Kjnt1f9+uvfTWC6hJc4hV3tAAoasK42+2Gy+VCW1sbjMYpbkCZhtPphCiKMJvNiX0GgwGiKCae+/1+NDY2oqOjI1HPYrHAYrGgs7MzPx+CiGgORKJRdO3/IpZVLcO9d//zrL9fWWk5ztSsRPSkd9bfi6jQFJU8mUwmmEwmtLe3Z32sw+FAZ2cnBEGYtI7NZkNdXV1KgmW322EwGOB0OlP2ExEVM/fff46bzz2HEYsdJRrNzIJUVQF33x0rMzC4/DosO304s1hZxCUqNoobtgOAurq6rOp7vV643W7YbDY4nU5IkjShjiRJcDqdMJlMKftFUYQgCHA4HLk0mYhozsiyjMFffBY9lUuwaSYTxeMaGoDHHouVGdCtM2LZ4EWMDfVOHSvLuETFRpHJU7b27t0LAGhvb4fFYkFtbe2E3iu32w0gNpQ3niiKideJiIrd3/66Fzedfx6aNz0EVYl25oEiEaC/P1ZmYPUNdwIAjh5+dOpYWcYlKjYLInmy2+2QZRk+nw8tLS0AAKvVmpJAxa+6SzesF+/pStdjNZloNIr+/v6URzAYnPmHICLKQDQSwYDzUzhdswpb7v5wbsGeeQaoqYmVGbjx+jvQry3DmWcfmzpWlnGJis2CSJ7iRFGE3W6Hx+OBIAiwWifeIDPdkGA8ocpmWYOuri7U1NSkPFpbW2fcdiKiTPxxfyvWBk6gcocd6pnOdZohbYkWF5ZdB1X33+f0fYnm2oJKnuKMRiPsdjuAqwlR/Kq7QGDiAm/xHqepJpuP19DQgL6+vpTH7t27c2s4EdEUhgd6UPbbL+Po2lfglbe/oyBt0Bhuw6rL3QgGRwry/kRzYUEmTwASE8PjiVE8eUo3NBdPqJKXNZiOWq1GdXV1ykOv1+fWaCKiKfzpex9CaSiIm973XwVrw7KbmlAaDeHF5zhPlOavBZs8xcXXi4onRj7fxHsz+f3+rNeVIiKaS+e7PVjldeKFLTtgWHdLwdqx4UYTRtUlOH/k8YK1gWi2KWqdp0z5/f5pe4nGr9skCALMZvOEq+r8fj8kSeKQGxEVLTkawYuO90BVWg3TP/5n/gLfdBNw8SKQxZSFstIKnF5UjxL/36aOlWVcomKiyJ6ndPOS4tra2hKLWgKxNZ4MBgPa2toSQ3JerxculwsdHR0px9rtdgQCgZQEymazwWQycYFMIipaT37/I1h6/gUE3vZF1FVmtw7elLRaYMmSWJmF4LotWHr2eUTDofSxZhiXqFgoKnlyu92wWq2QJAlerxdtbW3welNvBRBf1DLe8ySKIkRRRGtrK+rr62Gz2XDo0CG4XK4JE8BFUYTH44HD4YDNZoPFYoEoinC5XHP1EYmIsnLw4M+x9EkHnr3lbXhL08QriHPi8wHbt8fKLKzc8jZUhEdx9PAj6WPNMC5RsVDJsiwXuhHz0fbt27F///5CN4OI5rHL0nkc/eRNCJVW4zV7nodOV5rfN/B6gcZGwOMBspj3ORYew5EHanGu8T68+QM/mBgLmFFcormQyff3vJzzREQ038myjN//3x1oCA5iSYsr/4lTDnQlOpxdtQmlL/2p0E0hmhWKGrYjIqKYR5yfxib/nzG0/T9wTQGvrptM+ca7sbr3JC5cOlnophDlHZMnIiKFecGzH2sf3YOu616P27Y/VOjmpHXTq98FNYBn//j9QjeFKO+YPBERKYh0+giGv/UunKhdizs/9nOoVKrZe7OVK4GvfjVWZmnpimvRvageYe/+ibFyiEtUDDjniYhIIS6e64Kv9S4M6cqx/uO/QnlZ5ey+4bJlwMc/PuPDR268Bw1POjAycBll42PlEJeo0NjzRESkAGfO+/Hcl+4AoiGs+NivsGHlhtl/095eoLMzVs7ADXfugkaO4pDbkRorx7hEhcbkiYioyJ0578PRL92B6uAgVnzit9jY8Iq5eePubmDHjlg5A9euuxldS69D+M/fS42VY1yiQmPyRERUxM6cOYaXvvBq1AT7seQTj2BdwysL3aSsyK/5J6y87MfFF54EAISk8wVuEVHumDwRERWpM6eOoPuLr4U+HMTSFjfWXvfqQjcpa68xPYBT5Ysw8tPYHKdB7y8L3CKi3DF5IiIqQs95H8GZL74GMoBrdv8O6wxbC92kGaktq0TXPS0YVceuTzp2kUN1pHxMnoiIishQKIiOH/4r5IffhkB5HVZ98g9Yt/aWwjSmrAzYvDlW5qD5zR+H3/oDDCzV4dzQpbzFJSoULlVARFQkRsdG8eP/tOB1Rx/BhRuacOdH/xd6fUXhGnT99bF70uVIo1bjDW98B/5m/Qp6gkFEN1wHdR7iEhUKe56IiIpA8JIfT33qZrz26CMY2vYJvPrBRwubOM0CnbASS0YknB8eKHRTiHLC5ImIqMC6f/8dvPTJTVANXMLZf/o+Nr+zbXZXDs/U008Den2szEOsW2y/wcaTF3Dxj7/LX1yiAuCwHRFRgfRc9MPz7V1Y99Lv4Vl+E5a+5xt44/W3F7pZV8kyMDYWK/MQSxWOQB+N4tSZY/mLS1QATJ6IiOZYMBzCb3/+Baz7bRvqVBoc3vYgLPd9FhU6faGbNidGzh4rdBOIcsLkiYhojsiyDPfff4Hg3n/F9YET8Iu34ZZ//im2LF5d6KbNmb7yWmgu+grdDKKcMHkiIpoDR7v+isM/+gQ2n/w7equXIfrAT3HPrTsK3aw5N1y3FpW9pwvdDKKcMHkiIppFA5dO4E8/+CjWPvcIGnQVGHrzJ3HrWz4FtVYBQ3TXXw889xwginmLNfSn/4tFf/k+Qoefhvba63KPS1QATJ6IiGbByWNPwrf/C1j+whNYqtHBd/v70PR/voLS8upCNy1zZWXAxo15jVXVczvK//ptvKjtx41cJJMUiksVEBHlSTQaxZ9/9208YtuIoP1O6LoP4u+veBeELx7Bvfe3KytxAoCTJ4H774+VeYq1qW4jSgZlrLjv7vzEJSoAJk9ERDkKjg7hDz/7HJ782Gos+YEVurFhnHzbF7H5P8/gvR/4HgzL8jDsVQg9PcB3vhMr8xSrdEwFTRBY9EIQ8uXLucclKgAO2xERzUBUjuKg59c4/YfvYvWLv8OKsSG8uHITpHd8Fa+/dQfUav5tOpke08ew9tdfx6Xec1ha6MYQzQCTJyKiLFw4cRiH3N+E7tlHsLb/HNS6Clze8HoIb/xX3Lvh1YVuniKsXh7riTt57PdYanpzgVtDlD0mT0REUxgYG8XxF56E9LdO6J9/HIull7FarcXLaxsx9NZPo/G174W6RFvoZirK4vpXAABC3QcL3BKimWHyREQ0ztnec3A98V0MPvtbbLhwDKtHegGNDs9fcyPO32HF60zN2FS1uNDNnH3LlgEPPRQr8xhLBcC/pRojpUVw/z6iGWDyREQLXmSoFyNdf8HxZx5Fz9EnsOxiF26TIxisXIzghjsQNr4V1299OzaXVhS6qXNr5UqgtXVWYh1/zRIEKzgvjJSJyRMRLSjDoSCePvYnnH/eDe1JL+rOHcXSvrMAgFFtOS4uFhG684N45R3vRcOaTVCpFnDvyMAA4PEAjY1AVVX+YgGoPTOG82UjeWgk0dxj8kRE89pocBhHnnkULx9xQfb/HddcOo4lY0NYAuBc9XKcXXYduoxmqA2vwHpxK96xwrCwE6ZkXV3AnXfGkh6jMX+xAGzddxoH7lfYuldEVzB5IqJ5QZZlhPsu4Mzxp3D6xT9j6OTT0F94CUv7z6FGjqJUXYILSwwYMr4ddZvuhrjpHlxbWVfoZi9o0Wi40E0gmhEmT0SkGFE5iosDAVw4dQR9Lz+H4LkXIF/0obTnBKqlMygfGwYA1KhLMFi9Ar3LrsXoVjNW33AXrr/RhJt0Crif3AIiR5g8kTIxeSKiohKJRnDW78G5l59DT88ZDF7qhnz5BEqlM6gdvITFo/0og4wyAEMaHS5ULsEZYSWOr7wZ2hUbsHr9K3Hjhtuxqbym0B+FpiGz54kUiskTEc2Z8HAfLl3w4/JFH6Se0xjqPYtg71moe09D33cepSN9qBodQEUkCAGAAGC4RI++qmUI1q5EUHwFepc3oPqaG7Bk3c0wLBWhUWsK+6HmM602dpWcNg/rWI2LFazWI4JI7nGJCoDJExFlZDQcwsWRAVwY7seF4X4EBgIYHZYQHu1HZGQAkdEBYHQQuoGLqL3cDVU4CM3YCCpGJFSN9qFmdAClkTEAgB7AMgBhlRp9+ir0Vy7GaPUyjK66CSM1y1C5bgtWiFuwcmk9dKWVnMBdKDfdBLz88qzEOvyp12Hg7IsIRSPQMgEmhWHyRLTAyNEoosO9iIZGAVnG8GAvLl/oQv/Fbgz3nMJo/yX0asvQo9KgJHAapQMXEQ6NQo6EoIlGURcawprhAK6dZMglqlKhp3oFQroyREtKEaxdhUDlJlyqXAJ93UpULVqDuiX1WLpMxKLaFew5WqA0JXpo5CgGQ0HU6ssL3RyirDB5ypDf74coKvTO6FRwsiwjKssIy1GEoxGEo1GE5VgZiUYxFg4iMjqI8OggIqMDVx6DCF/ZjgaHEIpGMaavQrC0EqMlpYgGB4DhfmB0ABgdgHp0AJorD+1oPzShEUSgQhgqlIRHUT46gPLRflQEB6GRoxPaWAEgXFKGMW0Z1oeGcWM4iIGKRRitXgZNWRV0Oj102jLoa5ahctVGVNauhKa0EurSKqhKK6HWx7Y1VYuh1pXN/T8y5d+RI8Ab3gA8+mis5yhfsQBs+sIBVL12DJdPPYvahlvz0FiiuaO45EmSJLS2tsLr9cLlcmV8nNvths1mg9frhSiKsFqtaGlpSVvXYDDA7/en7PP5fFm189LIAB4/dTSrYyiJLEMVGoU6OAhVcBDq4BDUowOxMjgIVXAIUGsga3SQS7SQ1VqENSWIqEsQ0mgRUakxpi5BWF2CkFqD8JXnIbUaY9AgpNZgDEBElhGKRhCRo4iEQ1AHB6AZHYQmOIiS4BC0wQGUBIehCQehDo9CHQ5BEw5CEx5DSWQMcjSCKIAIVIhAjm3LAOQo1NEwVNEo1HIE6mgE+kgIZZExlEVCKI2GUJb0XH+lF0eF2A9ltj+YEZUaIyWlGNGWYkRXjhFdBUa0pVBDRokcwZiuAv3CSgTLBITKBYTKa6HWlaNGX4qqykUQltRj0TIDrhOWoUqrh0qlghyNQMVeoYUtFALOnImVeY6l6x3CgFbAyDffCdF+FBom3KQgikqe3G43XC4X2traYMxiwTan04nW1lbs3LkTW7Zswb59+2Cz2dDT0wO73T6hriiKMJvNiX0GgyHrXqfh7kM48f/M01fMkkqWoYIMtSxDBUAFGSpZhvpKqbpSR528f9y++HHqpFiADLWcuj/5fa7WiT2PqFSIqtSIQoWoSoWISo2oSnXleWw7gvj+K6+pVIggvn21vi4aRkU4iIrIGMojY7HtcBAayJP+OwTVJVDLUWjT9KBkKgoVwmoNIuoSqCCjNByctG5ErUFYo0O4RIeIRodoiQ6REh2g0lw5B7jy7xuFCoglduoSQBN/6CFrSyHrygF9BWRdOVT6CkBXgdHSCgT1FVd6bmI9OJrSKpSUVkJbXg1daTV05dXQAtCO9kM9LEEVHIa6rArqcgGashqo9BV5nxfExIlmW/hND6H2D/+On37lXrz9wUdQrtUVuklEGVHJsjz5N1SRqq2thSiK8FxZqXY6jY2NKXUlSUJ9fT0kScL4j9/U1ITOzk4IgpBTG003LEH7e27OKcZkVCoNoFYDUAEqdWxbFUuJoFZDpVLH9quu7E95jnGvj6ujVkMVq3Ql7pUHkp4DQDQSe8ixUo5GJ98nR4BIGPKV/UiqK8sRqLRlUJVVQ1VaNa6shqqs6moZf620CipN7IodORoFImNAeAyaaAiaSBjqSAiaaKyUw2OQw8FYGQpe2Q7GtiNjkEOx11UqNdQVAtRlNdCUC1BfeWjKBajLqqHSKOrvDKL88Hpjt1PJxwrjybGAxPavn/8Rrj3w/+G5ez+Nt9/3H7m3mShH27dvx/79+6eso8hvhLq6zFcF9nq9E3qXBEHAjh070N7enjKXyev1Job3mpqaYDKZZpxEla+/DaJt6n98IqKF7k3v/ioe6z4Iw2++hD+vuRm3b31LoZtENK15f0tro9EIk8k0YX88KUoejtu7dy8AoL29HRaLBbW1tWhvb5/R+0ajUfT396c8gsHJh4WIiIpOQwPwxBOxMp+xkrZVKhXuevBR9NQsh+a7/4TBnjwtjUA0i+Z98jQZr9ebMq8JAOx2O2RZhs/nS0wmt1qtM0qgurq6UFNTk/JobW3NS9uJiOZEVRXwutfFynzGGhdXW1aJxR/cB1UkhKdaX4/R0aHc349oFi3I5Mnv98Pv96OjoyPt66Iowm63w+PxQBAEWK3WrN+joaEBfX19KY/du3fn2nQiorlz5gywe3eszGesNHE3rX8F5H/6LpYHTsBtvxtylKuPU/FakMmT1WrNaFK40WhMzJcav3TBdNRqNaqrq1Meej1vSkpECnLhArBnT6zMZ6xJ4t56qxln3vZFGLr/ht9+7S0IhvOwRALRLFhwyZPNZoPFYsl4qYP4fClJkmaxVUREBABNb/4Ennn9R2B47lH88pvvxggTKCpCCyp5cjqdMBgMaG5uzvrYbNaVIiKimVGr1HjHP34NT9/6j7jF68T3/+u9hW4S0QTzMnlKN8TmdrsRCAQmJE5ut3vKWE6nc8LEciIiml13v/dhPLHhbrze8z/4ybfeh1H2QFERUWTyFAgEJn2tra0NBoMBTqczsc/tdsNqtUKSJLS1tSUeVqsVDocDQOzqO4PBgLa2tsQQXfwWMJNNLCcimtcWLQLe//5Ymc9YGcStLa3Arpbf4JnN92HLUz/AH37wkdzbQJQnilok0+12o7OzE5Ikwev1oq2tDSaTKWVITRRFCIKQsvBlU1MTgNh8p/E6OzsTx4miiNbWVrS2tqK5uRkGgyGr++cREc0ra9cC3/727MTKIK5arYb5I3vxk/+7E1ufbMcjwRG84YH/zvutiIiypcjbsyhBJsu7ExEVtZERwO8HRBEoy/HGvcmxgKziDofG8KNv/gNe97QTT4m3w2z7Lcr15bm1h2gSmXx/K3LYjoiI5sCxY8CNN8bKfMbKMm65Vof3ffgneOK1D2Cr/y/41X+8Eie4EjkVEJMnIiIqelq1BtZ/+gaC7/sONlzqwvOfvQ2+k88Wulm0QDF5IiIixbjljvdA9ZGfQwgO4PIXX4Nf/e67iESjhW4WLTBMnoiISFE23fIGrPyPv+J8eR3qf/gAfr33U4VuEi0wTJ6IiCg9lQrQ6WJlPmPlIe66lddjm/15vHTtHbj+sTY8/G+vwMW+i7m3kygDilqqgIiI5tDmzUAwODux8hC3VF+Otz70GPY8/B68+el9OPKpm6F7XwdebXwTlzOgWcWeJyIiUiy1So1PfviHUH38UURkGcI37sO+H3wMkWik0E2jeYzJExERpXfsGGA05m+pgnisfMa94qYbX49Xtb2AZ8VXYfMTD+NX/74VL548krf4RMmYPBERUXojI8DTT8fKfMbKZ9wklRUC3vVvT+CEuQ0rLx3HwOdvxZ5v3Y9QJJzX9yFi8kRERPPK3W/+BNZ88Vn8ra4eb3/qv/HrT23G011/K3SzaB5h8kRERPPOsiXr8IZPPQG/5ctY0XsK6tY78InP3okzfZcK3TSaB5g8ERHRvCTWLME9b/o4Vn/hCB5ZsQn3d/8Rxx66AQ93fhYj4bFCN48UjMkTERGlV18P7NsXK/MZK59xM7ByyRq866Hf4uf3fRmDpVW4+zefw08e2oT9f//FnLw/zT8qWZblQjdiPsrkrsxERDS3ZFnGv3/zfbjn2V+gbmwIh8TbEb37Y3jPK97CtaEIQGbf3+x5IiKi9C5cAL72tViZz1j5jJsllUqFz3zgOzj/0d/Atekt2Hjy79jseAf+7TOvwV+6n5nz9pAyMXkiIqL0zpwBPvGJWJnPWPmMOwMlag3MN9yOj/6LE+c//jh+sXIz3nnq79B/8TZ86tOvxifd3+PNhmlKTJ6IiGjBetMNr8anP/9XPPG+7+GxZTdi5+mD2PHTB/CFT9+OLz7xAyZRlBaTJyIiWvA+/Nr/A9vnn4Lvo7/BL64x4i1nn8a9P9yFr/7bVnz/Dz/E5ZHBQjeRigiTJyIiIgAatRpvu9mET3/uzxja/Sf8z+pX4LWXXsJt//1ePPqpW/C+h/8Jv/I/g6FQnm6WTIrF5ImIiNKrqQHuvTdW5jNWPuPOApVKhdvXb8GXPvtHlH/+MPZc9wbUjQ1h96HvQ/+VJnz94Xdj//NPFrqZVEBcqmCWcKkCIqL5YTAUhF+6CNcfvgfhqR/jVZePI6zW4MnF16LkVf+AV7zqHdi8dE2hm0l5wqUKiIho5kIh4NKlWJnPWPmMOwcqtXpsWrIanzD/O3Z86TAu2X6PAxvfhOsHzuLOX+xG76c348HP3oknn/s9QtFIoZtLc4DJExERpXfkCLB0aazMZ6x8xp1jVbpSvPb6V+Njn/gZGv/zZbje8kUcrb4G/3DyL1j+lbvwy38R8YX/3AnfmRcK3VSaRSWFbgAREZES1ejL8MG3PYThN38cI0MB7Pnev2DTqYO47/D/Yuzw/+KHtevw+6UbcPc9H8G2G14NQV9e6CZTnjB5IiIiykG5VodyYTm+/LGfIhKN4sKlE/j5zz6PZS/8Dg+++CjkF38LV81K/HnRehxZcSN+/K7P49CFE2hcuhYrKopz0jxNjckTERFRnmjUalyzTMQHP/Df6BkdxNETz+BHP/sCXnX5OO7v/iN0/t/j8OH/wcE6A75ZV4897/w8Ni0XC91syhKTJyIiolmwqLQSr9lwO2576Nc4P9yPj7j+G/quP+H2ni686fyzeNfpv2HkyP/i+8IaHKpbh0O16/Chu96PbWs3IhSNQKvWYHFZZaE/BqXB5ImIiNK7+Wagrw+oqMh/rHzFVYAStQarKmvxs7d9HMDHERgdws0/+RzWD17Elt4T2NJ7Ag/4noBWjuLCs534ibAGz9aswrPCKnRY/g03L1ld6I9A4zB5IiKi9DQaoLp6dmLlK64C1ZVW4PQ/2RPPo3IU7/n1NzD84h+xpfcENvWdxt0XnoMaQM/hn8JRswpHalZB33A7xOtejVuvWY/GpWsL9wGIyRMREU2iqwv40IeAhx8GGhryFwvIX9x5QK1S44f3fhi498OIRKM4fPk03vLLr+PGvrPY1HcaN/W9jA/4noD2+AGMPKbFS1XL8EjVCrxQtRy1DbfjC2/8IKp0pdCoufrQXGHyRERE6Q0MAI8/HivzHStfcecZjVqNxqVrcfT9X0NYjkKvKcHvXn4Rvzvvw6N//DGuHziH6wbO445LL2LnyweBY7/Ci49+Di9ULcdLVcvgr1iC4NJr0b7jU1hcvnB792YbkyciIqIio1GrobmyjvXrV12H16+6Dh/e8kZE5ShC0ShC0Qj+9+nH8Ksnf4gN/eewYeAc7j37DOpCw8BR4MyT/x/+VLEY/sol8FfEHtde/1rU1q3EP264FfXVi1PeT5ZlqFSqQnxURWLyREREpBBqlRp6jRp6TQne84rteM8rtgOI3X/vu0f/jKMnnkXlZT8iZ56HOHQJG/rPYdv556CVo8Az/wNJWwZPWR1+Ub4Ip8trcap8EU6V1+F8aQ1KtaX41Zs/iP3dz6BGV4ZdG18NABiLRqDXMF1Ixn8NIiIihavU6vGRm18P3Pz6lP0Xhwfw02N/wZ8PPwbtZT/WDAeweiSAhsELeP3FYyiLxu4vGFKpcaasFk89/VNES2twtEzAOx6twbnSGpwrFTBSosPtKwzQqNR48mwXfv7GB7B12ToAwKWRASwpq5rrj1xQTJ6IiCi91atjk7pX5+FS+fGx8hWXprS0vAofbdyGjzZuS9n/9wsn8KE/dsIghxA4/SxWDweweqQXK0b6sFk6hTecP4LSaDhRv6+kDOfKYsnUDfoq/Kj9EL6ur8IlfRUu6yrRo69EQ+0KvHJ5PX7wwlO4ZfFqtN3+dlxfuxwqlQrnh/uxfJI5WH3BEVRq9Yqa8K6SZVkudCPmo+3bt2P//v2FbgYREdG0onIUP3zhb1hcVokfvPAU/nz2OGpDw1gxImH5aB9WJD0WBwewJDiQklxFoEKvrgKX9ZW4dCWpCmgr0KsrR6+uApK2PLatLceYRpu2DQfe+i+o0upxrPc8uqSL2F6/CWE5iiqtHrWlc7cmWCbf34rreZIkCa2trfB6vXC5XBkf5/f7YbPZIIoiJEmCIAiw2+0zrkdENO8FAsAjjwBvfCNQV5e/WED+4lJeqFVqvOf62wAAb1p3U2J/l3QRv+x+Bh/edCd+8uLf0fL3XyMiRwFZRlV4FIuDg1gSHMCSsQEsDg4knht7T6J2bBg14ZEJ7zWk0V1NpnQV6NWWo09bjv96+BD6S0oxUFKKfm0pvq8tw8CV5xG1JiXGvnt24VUrDLP7jzIFRfU8ud1uuFwutLW1wWg0wuPxZHSc3+9HY2MjOjo6YDabAQAWiwUA0NnZmXW9TLDniYgUz+sFGhsBjwcwGvMXC8hfXJpzUnAYgr4cABCJRvG1w27880134Lmes+gfG8H7DvwgUVcTjaAmNILa0BBqx4ZRGxpG7dgQhKTt2tAwqkMjqA6NJuZgjTek0cUSKW0p+ktKcaz6Gnzuc3+elc8373qeTCYTTCYT2tvbszrOZrOhrq4ukRABgN1uh8FggNPpTOzPtB4REdFCFU+cgNiSCg8a7wYAvHJ5PQDg5fftSakvyzK6+3sQkaO4NDKALUvXwnvpFP7fM0/gD2e7Uupqo2FUhkdRHRpFVVJZFRpBdTi+PTrLn3B6ikqe4uqy6OaVJAlOpxPNzc0p+0VRhCAIcDgcMJvNGdcjIiKizKlUKog1sXWlGoSlAIBbl4u4dbk46TFScBiei6cwHB7D1mXrcPDCCZwZlHBioAff9R+G8w3WOWn7ZBSZPGXD7XYDAAyGiWOjoigmXs+0Xqai0Sj6+/tT9un1euj1+qziEBERLTSCvhx3rd6QeH5v/abE9p5Xva0QTUqhnOsCZ8jv9wMABEGY8Fq8B0uSpIzrZaqrqws1NTUpj9bW1uwaT0RUSBUVwK23xsp8xspnXKICmPc9T3HphvriiVI8ccqknjHDyY0NDQ04ePBgyj72OhGRolx3HfDXv85OrHzFJSqAeZ88iWJsTDUQCEx4Ld6TJAhCxvUypVarUV3NmzISERHNN/N+2C6eFKUbcosnSqIoZlyPiGjB8HoBlSpW5jNWPuMSFcCCSZ58Pt+E15KH4TKtR0RERAvbvEyekucwCYIAs9k84Wo5v98PSZKwe/furOoRERHRwqbI5CndvKS4tra2xKKWcXa7HYFAICUxstlsMJlMExbEzKQeERERLVyKSp7cbjesViskSYLX60VbWxu848bM44taJs9PEkURHo8HDocDNpsNFosFoihOuDdepvVmyze+8Y2iiJfpcZnUm6pOtq/l+98nX3jepn6tWM8bkN+28bzNnYV03qZ6Pdv9haak8zYtmWbFvffem/Ux119/fV7bMNN4mR6XSb2p6mT7Wr7/ffKF523q14r1vMlyfts2L8/byIgsd3XFylwlx8ox7kI6b1O9nu3+QlPKecvk+3veL1VAREQzVFoKrF8/O7HyFZeoABQ1bEeULBgM4jOf+QyCwWChm0JZ4HlTkO5u4N3vBrq7cz9vSbFStmnW8Wcu/5g8kWIFg0F89rOf5S8EheF5U5DeXuDHPwZ6e3M/b0mxUrZp1vFnLv+YPBERERFlgckTERERURaYPBERERFlQSXLslzoRsxHGzduhMFgyOqY7u5u1NfX560NM42X6XGZ1JuqTravjd8XjUbR1dWFhoYGqNWF+zuA523q14r1vAH5PXc8b3NnIZ23qV7PdH+xnDulnDefz4fnn39+yjpMnoiIiIiywGE7IiIioiwweSIiIiLKApMnIiIioiwweaIFwev1orGxESqVCk1NTZAkqdBNoiw4nU7U1tZOuBE4FQdJkmC1WtHe3g6r1Qq/31/oJlGG+LM1M0yeaN6TJAmtra3o6OiAy+WC3++HxWIpdLMoQ5IkwWQyMeEtYhaLBRaLBc3NzbDZbGhqaip0kygD/NmaOd4YmOY9t9uNjo4OCIIAALDb7UyeFCR+3qg4eb1euN1uuFwuAIAoigBiPRpms7mQTaNp8Gdr5tjzRPOe2WxO+SUhCAKMRmPhGkQ0j7jd7gk/T6IoJpIpovmIPU9UMPHhNK/XO+kvWr/fD5vNBlEUIUkSBEGA3W7P6X1dLhd27tyZU4yFrlDnjvIrH+fR5/Ohrq4u5RhBEDjvaRbx568IyEQF4HK55JaWFhmAbDQa09bx+XyyIAhyZ2dnYp/ZbJbNZvOM37e3tzen46lw5w6A7PF4Znw8pcrXeWxubpZNJlPKcen2UX7Mxs8ff7ayx+SJCkoQhEl/AZjNZlkUxZR9Pp9PBpD4pdDb2ysbjcZJHz6fL+X4lpYWube3d1Y+y0Iz1+eOv+BnR67n0W63TzjebDbLLS0ts9NgkmU59/OWjD9b2WPyRAUlimLaXwC9vb0yALm5uXnCa4IgzOivWofDMeELmWZuLs+dLPMX/GzJ9Tx6PB5ZEIQJMdN9SVP+5PPnjz9b2eOEcSpKbrcbANLeXFkUxcTrmXI6nYl5GX6/H16vF+3t7bk3lCbI97lLFggEZnwsZSfT82g0GiGKYmKdIL/fD0EQeKVdgcz0548/W9nhhHEqSvHJpukupY0nQfFJkNNxu91plybw+Xw5tZHSy+e5i9eNJ7qdnZ0QRTFxOTzNnmzO44EDB2Cz2WAwGODz+dDZ2TmXTaUk2Zw3/mzNHJMnKmrjr+IBrv5S8Pv9GS05YDKZIMtyvptG08jHuYsf09LSgpaWlnw2jzKUyXkUBAEOh2OOW0ZTyfS88WdrZjhsR0Up/tdPuq7k+Gq4XOCtOPHczQ88j8rE8zY3mDxRUYr/Akh324D4LwV2Lxcnnrv5gedRmXje5gaTJypK8R/udPOSshnyobnHczc/8DwqE8/b3GDyREUpfrXO+CtD/H4/JEnC7t27C9Qymg7P3fzA86hMPG9zg8kTFdRUl8fa7XYEAoGUXwI2mw0mk4mXQRcBnrv5gedRmXjeCksl8zIkKgC3243Ozs7EZbJ2ux0mk2lCl3Ly/Zn8fj9EUeT9mQqM525+4HlUJp634sDkiYiIiCgLHLYjIiIiygKTJyIiIqIsMHkiIiIiygKTJyIiIqIsMHkiIiIiygKTJyIiIqIsMHkiIiIiygKTJyIiIqIsMHkiIiIiygKTJyIiIqIsMHkiIiIiygKTJyIiIqIsMHkiIpoBq9Va6CYQUYGUFLoBRERKY7FY4Pf7C90MIioQ9jwREWXI7/ejsbERTqcTJpOp0M0hogJh8kRElAG32w2LxQKv1wsAcDqdaGxsTDwnooVDJcuyXOhGEBEpRW1tLSRJAn91Ei1c7HkiIsqQ3++HJEkcsiNa4Jg8ERFlyO12AwCampoK3BIiKiQmT0REGXK5XAAAs9lc4JYQUSFxzhMRUYZqa2sBAL29vQVuCREVEnueiIgy4PV6Od+JiAAweSIiygjnOxFRHJMnIqIMxOc7mUwmuN1uru9EtIAxeSIiykC858ntdkOSJBiNxgK3iIgKhckTEVEGTCYTBEGAJEm82o5ogePVdkRERERZYM8TERERURaYPBERERFlgckTERERURaYPBERERFlgckTERERURaYPBERERFlgckTERERURaYPBERERFlgckTERERURaYPBERERFlgckTERERURaYPBERERFlgckTERERURb+fyHWNtgK7udnAAAAAElFTkSuQmCC",
+      "text/plain": [
+       "<Figure size 640x480 with 1 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "plt.plot(t_hist, Ft_hist, label='histogram')\n",
+    "plt.plot(t_cont, Ft_cont, label='contour')\n",
+    "\n",
+    "for imtd in lens.img_td:\n",
+    "    plt.axvline(imtd, ls='--', c='r')\n",
+    "plt.legend()\n",
+    "plt.xlabel(r'$t$')\n",
+    "plt.ylabel(r'$\\widetilde{F}(t)$')\n",
+    "plt.xscale('log')\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5df8eb60-c0f6-48a3-b8d5-0099c492f224",
+   "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.10.14"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}