Skip to content

Uddeepta Deka / Updated Charged Lens

Go to a project
Commit 02a536ab by Uddeepta Deka
Options

python script to generate Fw

parent 572f5382
Inline Side-by-side
Showing with 384 additions and 331 deletions
notebooks/Ft_computation.ipynb deleted 100644 → 0
{
"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",
"- [Approximate Ft](#Ft_approx)\n",
"- [Testing $\\widetilde{F}(t)$](#Ft_test)\n",
"\n",
"This is implemented in [generate_Ft.py](../utils/generate_Ft.py)"
]
},
{
"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 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",
" # We assume that atleast one image exists\n",
" assert len(lens.img)!=0\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",
" # We assume that atleast one image exists\n",
" assert len(lens.img)!=0\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",
" if len(level[i])!=0:\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": "732fb49f-c634-4228-80f3-2ea22231eb57",
"metadata": {},
"source": [
"## Approximate $\\widetilde{F}(t)$ <a class=\"anchor\" id=\"Ft_approx\"></a>"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "e819b1b6-c405-46cc-87c0-feee4e4f4ffe",
"metadata": {},
"outputs": [],
"source": [
"def approx_Ft(lens, t_max, dt):\n",
" \"\"\" \n",
" Computing F(t) due to the stationary points.\n",
" Ref: Eq. 5.5, Nakamura, Deguchi, 1999\n",
" \"\"\"\n",
" assert len(lens.img)!=0\n",
" image_timedelays = lens.img_td - min(lens.img_td)\n",
" t_arr = np.arange(0, t_max+dt, dt)\n",
" Ft = np.zeros_like(t_arr)\n",
" for ii in range(len(lens.img_td)):\n",
" img_morse_idx = lens.img_midx[ii]\n",
" if img_morse_idx == 0.5:\n",
" Ft -= (np.sqrt(np.abs(lens.img_mag[ii])) * \n",
" np.log(np.abs(image_timedelays[ii]-t_arr)) / np.pi)\n",
" elif (img_morse_idx==0) or (img_morse_idx==1):\n",
" sel_idxs = t_arr > image_timedelays[ii]\n",
" Ft[sel_idxs] += (-1)**(img_morse_idx) * np.sqrt(lens.img_mag[ii])\n",
" return t_arr, np.asarray(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": 7,
"id": "9a43cfbe-15a6-482b-a611-b3af5571132a",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Tmax = 28.835924899983453 | dt = 0.020943951023931952 | xmax = 9.113038116454486 | dx1 = 0.0034570190957158614 | dx2 = 0.002312106639417798\n",
"Number of td samples = 41567051\n",
"x_max = 9.113038116454486 | t_max = 28.835924899983453 | dx1 = 0.025, dx2 = 0.025\n",
"Number of td samples = 1432\n",
"\n",
"Time taken to compute Ft_hist 6.553 s. and time taken to compute Ft_cont 10.116 s\n"
]
}
],
"source": [
"qtest = 0.05\n",
"ytest = 3\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",
"t_approx, Ft_approx = approx_Ft(lens, T_MAX, dt)\n",
"\n",
"print(\"\\nTime 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": 9,
"id": "21f428d3-2346-4d5c-b94c-0fd503e51900",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAkUAAAG6CAYAAAAGUjKQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAACB+ElEQVR4nO3dd3hT1+H/8bcs7ylP9rCM2QSwIZA9sJOmbUjS2NCk6UgacPf8Fof21zZdIaZ715A0bZO0BbtpkzbTIpsssNkby2wwYFveW/r9cS3ZwgaMGZLw5/U8fu710ZV0rHOxP5x77jkml8vlQkRERGSQC/J1BURERET8gUKRiIiICApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAgAwb6uQKCZMmUKaWlpF/ZFW1pg/34YMwbCwy/sa8ul4W7DoUPh2DGjLUHtKiLiJ8rLy9m2bdsZj1EoOkdpaWk899xzF/ZFy8ogMxP+8Q/IyLiwry2XhrsNv/c9uO8+oy1B7Soi4ifmz59/1mMUivzBiBHw858bWwlM7jacPt27LdWuIiIBQ6HIHwwZAt/4hq9rIeejZxtOndpdrnYVEQkYGmjtD2pqoKjI2EpgcrdhRUV3W6pdRUQCikKRP6iogAULjK0EJncbvvNOd1uqXUVEAopCkYiIiAgKRSIiIiKAQpGIiIgIoLvP/ENEBMycaWwlMLnbMC7Ouy3VriIiAUOhyB9MmmRM/ieBq2cbfvSj3eVqVxGRgKHLZyIiIiIoFPmHDRsgLMzYSmByt+E//tHdlmpXEZGAostn/sDlgrY2YyuByd2GTqd3W6pdRUQChnqKRERERFAoEhEREQEUikREREQAjSnyD5MmwdatYLX6uiYyUO42HDbMuy3VriIiAUOhyB9ERMCUKb6uhZyPnm2YkNBdrnYVEQkYunzmD/bvhwcfNLYSmNxtuHZtd1uqXUXkFMXFxWRnZ2Mymc56bFlZGfHx8ZRpEthLRqHIH1RVweOPG1sJTO423Levuy3VriJyipycnH4fW11dPaD3sNvtA3qe6PKZiIjIJWXt5zjDrKwsampqzvn1s7OzKS8vP+fniXqKRERELhvZ2dnqKToP6ikSERG/0dzRxl7HCV9X44zGWZKJCA4979dxOBzk5+ezevVqAFauXOl1ec3hcLB69WqsVitZWVmAcWksPz+f2bNns27dOmw2m6c3qbi42BOI8vLysFgsFBQUAGCz2SgqKsJisVBWVobFYmHlypVYLBavOrmPA1i9ejUOhwOLxYLVamXNmjUArFixgsLCQkpLS1m0aBE2m401a9ZgtVrJz8/HYrHgcDiw2+0UFBSQkZHhqXtxcTGFhYUUFRV5fpbq6mqWLl3KkiVLKC4uZtmyZZSVlbF48WIKCwvP+3M+FwpF/mDIEHjoIWMrgcndhhMmeLel2lXknOx1nOC2//7W19U4oxdv/zLTkkac9+vk5+eTm5tLfn4+2dnZLFq0yBOKysrKyM/Px2azUVBQ4AlFubm5LF261HNcdna25/VycnJYtWoVdrvdK0ysWLGCgoICr0tq2dnZpKamUlFR4QlGZWVlnktvVquVvLw8MjMzycrK8gQlm81GYWGhVziz2+1UV1d7Ak5paSlgBLN58+Z5QltCQgLl5eXY7XaWLVtGXl4epaWl5Ofnk5+fT3l5OdnZ2axZs4bVq1eTl5dHdnb2OY3DOl8KRf5gxAhYtszXtZDz0bMNZ83qLle7ipyTcZZkXrz9y76uxhmNsyRfkNcpKCjwBJKcnByWL19OWVkZGRkZZGRkUFBQQGZmptdzysrKWLdunScouHuCTsfdG7V06VKv8sLCQtLS0li0aJEn8LiDlHvMU0ZGBlar1etyXFZWlqeuubm5ZGVlsWTJEs/z3b1CAGlpaTgcDk9vk8ViITs7mxUrVpCXl+cV9FasWIHFYvH8XAsWLCAvL4+SkhKFokGnvh5KSyEzE2JifF0bGQh3G44fD7t3G20JaleRcxQRHHpBemECQc9LV7Nnzwa87zg79dIWdIenUy9Nnc769etxOBy9BndbrVasVis2m63X+7lDDBjByOFw9Pnas3r+BxA84QqMS3klJSVnrFvPugAkJib2qstA78AbKA209gd79sBNNxlbCUzuNnztte62VLuKyAVWVFTkGXuTlpbGihUrzni8u5enr3BhtVq9Ao97HFJ+fr6nzH35rr/cvUBWq5Xc3Nx+P89fqKdIREQkgBQUFJCdnU1eXh55eXkkJCSc9hKTuxfmdLfo9+yNslqtlJaWesY5AZSWlvZ7CoHs7GwsFounx2j9+vX9/ZH8hnqKREREAsTy5csBY2yPe0D0qlWrTnt8VlYWFoulzx4lu93O4sWLPd+7xx+VlpZSUFBAQUFBvwNRWVkZNpvNcxnQ/Xpw6S+BnQ+FIhERkUvIHRJ6DmB27/e8nOXer+oxK/6pt6hbLBavIOIOMcXFxZSVlWG321m5ciUOh8MrGBUXFwPeA7VXrFiBzWYjPz+f5cuXe756jjs6tY496+GuX1lZGStWrPCMKbLZbJ6lSs4UkPrqzTrdeKaLxiXn5Pbbb7/wL7p5s8s1YoSxlcDkbsPi4u62VLuKyCmKiopcFovFBbhycnJc5eXlrpKSEpfVanUBroyMDE9ZVlaWC3BZrVZXYWGhy+Vyeb4vKChwLV682LVkyRKv1y8vL3dZrVaXxWLxPMflcrlKS0tdWVlZrpycHNeSJUt6Pc/lcrlqampcGRkZroyMDE8d3V/u4wsLC73qX1JS4nn+kiVLPD9DUVGRq6amxmW1Wr1+poyMDBfgysrKcpWWlrrKy8tdOTk5LsBlsVg8z1u8eLFX2YXQn7/fJpfL5bq0MSywzZ8/n+eee87X1RAREbmgli9fTkZGhudWeTe73e65PT6Q9efvtwZai4iIDHLuySL76iexWq1ek0RezjSmyB9s2QIjRxpbCUzuNvzXv7rbUu0qIgHCPdYnNzfXs1yIe1mOvLw8rwHZlzOFIn/Q3g6HDxtbCUzuNmxp6W5LtauIBAj3Uh52u53c3FwyMzPJy8ujurqawsLCPieSvBzp8pmIiIiQk5NzSZfU8EfqKRIRERFBoUhEREQE0OUz/5CebqyZlZ7u65rIQLnbcPx477ZUu4qIBAyFIn8QEwM33ujrWsj56NmGw4d3l6tdRUQChi6f+YPDh2HpUmMrgcndhuvXd7el2lVEJKAoFPmDykp49FFjK4HJ3Ya7dnW3pdpVRCSgKBSJiIiIoFAkIiIiAigUiYiIiAC6+8w/JCbCZz9rbCUwudtw7FjvtlS7iogEDIUifzBmDDz2mK9rIeejZxtec013udpVRCRg6PKZP2huhm3bjK0EJncbVld3t6XaVUT8iN1u93UV/J5CkT/YsQOmTjW2Epi62vC5PxR0t6XaVUT8SHZ2tq+r4PcUikQuoFcObvd1FUREesnOzlZPUT8E3Jgih8PBsmXLKCsro6SkpN/Ps9ls5OfnU1ZWhtVqJS8vjyVLllzEmoqIyLlytjbRdnSnr6txRqHDJhIUFnler1FWVkZhYSEA69evJysri4KCAq9jbDYbRUVFWCwWysrKsFgsrFy5EovFAhiXwwoLCykuLqa0tJRly5axYsUKEhISKCwsJCsrC4Di4mJPIMrLy8NisXje62zvsWLFCvLz83E4HLhcLgCWL1/OsmXLcDgc1NTUeI4rLCyktLSURYsWYbPZWLNmDRkZGef1OV1qARWKbDYbJSUlLF++/Jw+6OLiYpYtW8bChQuZNWsWq1evJj8/n6qqql4noYiI+E7b0Z0ceHi2r6txRqMfXkf42IH/sbfZbBQUFHj+Y2+z2TyXttx/k1asWEFBQQHl5eWe52VnZ5OamkpFRQUWi4WEhATKysqw2+0sWrSIhQsXsnDhQnJzc8nNzfUElpycHFatWuUJUW79eY/FixdTUlJCcXGx55glS5ZQXl7OihUrACPUFRYWYrfbyc/PZ/bs2djtdqqrqwf8GflKQIWirKwssrKyPA3RX8uWLaO0tNTzfUFBAampqSxfvtw/QpHJBKGhxlYCU1cbOjHRGmwmzN2WaleRcxI6bCKjH17n62qcUeiwiQN+rsPhIDc3lzVr1njKsrKyWLx4MbNnz/Yck5+fz9KlS72eW1hYSFpaGosWLfL07mRnZ2Oz2Vi6dKmnsyAnJ4fly5djt9uxWq2nrUd/3gMgISGh1/PdPUnu+rvfMzc3l6ysrIC9EhNQocitrwY6nbKysl7Bx2KxsGDBAlasWHHGk+aSmTkTWlt9Wwc5P11t+NwTD/Hc3MkcmjnTKFe7ipyToLDI8+qF8Xfr16/H4XD0utrRswfHfcypf5usVitWqxWbzdbrdXuGlLS0NMAIPmerx7m8R3/MmjVrQM/zFwEZis7F6S6zuU+gcw1ETqeTuro6r7KwsDDCwsIGVD8RERk8ysrKznqMe/xPX5efrFbrBRkwfSneIxAN2rvPysrKyMnJOefn7dmzh7i4OK+vZcuWnV9lduyAjAzduh3Iutrw5o17efHhv3Tfkq92FZEe3P8hP1NPjPs/6z3H+vT1GufjUrxHIBqUochut2O321m5cuU5Pzc9PZ3a2lqvr1OvyZ6z5mbYsEGT/AWyrjaMbWll2oHK7skb1a4i0oP7jrD8/HyvcvcYH/cxFoulz/GzdrudxYsXX5B6nOt79Lwcd6ZLc4FsUIaivLw8zyC1cxUUFERsbKzXly6didvIxsC720JELh2r1UpOTg5lZWWkpaWRn59PXl4e8+bN8/oP9sqVK3E4HF6hxX0HWM9xslVVVYB3SOmr98fdM1RcXOy5Y62/75GZmQngudU+Pz/fc3nN/dzLJSQNulCUn59Pbm5uwM2dIIHh0wfe8XUVRMTPFRUVee7OWrFiBdXV1b3+o56Tk0NpaSlFRUXk5uaSn5/PunXrvAKPzWbzhJhly5Zht9u9ynqGl7y8PKxWK4sWLWL9+vWecHa29wBYvHgxOTk5FBcXk5eXR3Z2Nrm5ueTk5GC1WlmxYgWrV68GuoNToDK53LMxBZC0tDQsFovXbfb9UVxcTHV19Xl1Pc6fP5/nnntuwM/vU1kZZGZCaakxBkUCT1cbHr0Ohr2F0ZagdhUR8RP9+ft9WfYU9TVq3maz9RmI/CLRpqbC6tXGVgJTVxs2J8ORG0zG92pXEZGAEpC35FdXV592PNDy5cvJz8+nqKjIc3eZzWYjLy+PvLw8li9f7jm2vLyc6upqz8A3n4mPh9xc39ZBzk9XG3Y8H0RDTNf3oHYVEQkgARWK3Gu0OBwOysrKWL58OVlZWV7jg6xWKxaLxTOorKyszDN9+qmj/QHPjJ0+VVkJTz8Nn/gEDBni69rIQHS1YWiNk8gjJuN7ULuKiASQgBxT5EsaUyR90pgiERG/NmjHFImIiIicK4UiERERERSKRC4KXZUWEQk8CkX+IC4Obr/d2Epg6mrDjnBoGAnO2Fi1q4hIgAmou88uW2lpcKEHb8ul1dWGzZ8x0zwcUtOsmIPMalcRkQCiniJ/0N4OJ04YWwlM7jbscGJuceFsa1O7iogEGIUif7BlC6SkGFsJTF1tGHMA0la5cG3eonYVEQkwCkUiF4ELDbQWEQk0CkUiF4HT5fR1FURE5BwpFIlcBM7OTl9XQUREzpFCkchF4HRqcLWISKDRLfn+YPp0qK2FqChf10QGqqsN678QR+M9JiyTxoNlqNpVRCSAKBT5A7MZYmN9XQs5H+42DA7CCTidnWpXEZEAo8tn/mDPHrj1VmMrgamrDSOOOBlR4oTdu9SuIgPQ1tbJoYO1fv3V1nZ+YwYdDgd5eXnk5+eTl5dHdnY2ZWVlXo8vX76ctLQ0z7Hx8fGkpaWxYsWK8z4uNzeX+Ph4z3vabDZPfbKzs8nNzcXhcHieX1ZWRnZ2NiaTibS0NOx2OzabDZPJRF5eHna7/bw+D3+iniJ/UF8Pr7xibCUwdbVh8HUQdQQa6xwQoXYVOVfHKxv4zS/e8nU1zugr37iOkaMGvnxPbm4u1dXVlJaWApCXl8e8efOoqakBYP369RQWFmK3272CSl5eHnl5eSQkJJCTkzPg42bPno3dbqe6upoVK1ZQUFBAeXm5p37Z2dmkpqZSUVGBxWIhIyODkpISTyBKSEigrKyMxYsXU1hYeH4fpp9RKBK5CFydHb6ugkhAShkSzVe+cZ2vq3FGKUOiz+v57qDh5u7BcTgcWCwWsrKyyMnJYfny5RQUFGCxWAAoLS0lPj6e/Px8cnJyBnRcbm4uWVlZLFmyxNNrtHTpUq/6FRYWkpaWxqJFiygqKvKUFxUVkZmZSX5+PuvXr/eEusuJQpHIReDs1N1nIgMRGmo+r16YQNAzaBQXF1NSUtKv57kDk81mG/Bxs2bN8uyvX78eh8OB1Wr1OsZqtWK1Wns9PyMjg8WLF7NixQqvn+FyojFFIheBeopE5ExWrFhBXl4eVquV3Nzcfj/v1ABzPse5xwJVV1f3+fye44rc8vLyAFi2bFm/6hFoFIr8wahR8LvfGVsJTF1t2JIAlXNMtCUnqF1FpE/Z2dmUlJRQWFjodRmtP6qrq/sVePpznPvxnuOJenJfjuspPz+fJUuWUFZW5jWY+3KhUOQPkpPhi180thKYutqw3RJE7UQTnXERalcR6aWsrAybzcbs2bM9Ze4emb56bPp6fk5OzgU5LisrC4vF0me4sdvtLF682KvMPSapoKCAjIwM8vLy+uxNCmQKRf6guhqeesrYSmDqasPgOicx5S6oqla7ikgv7t6XwsJCT2+Le0yRzWbzujXfXebmDi8FBQW9Xvdsx50uvKxcuRKHw+EVjIqLi3s93263U1JS4glKK1euBDinS3+BQKHIH+zbB5/8pLGVwNTVhhEnYdjbLsyHDqtdRaQXq9XKkiVLsNvtLFq0iISEBIqKirBarRQWFva6ZLVu3TrPrfalpaWnvePrTMetWLGC1atXA7Bo0SKvAJWTk0NpaSlFRUXk5uaSn5/PunXrvC6p2Ww2srOz+7ycZrPZes1rFMh095nIReDs0N1nItK3goKCXr09pxvXs3Tp0j7DyLkct3jx4l6Xwnpyz0N0OllZWb3ql5GRgcvlOmu9Ao16ikQugvaOVl9XQUREzpFCkchF0NnR5usqiEgA6+/lqMvlspW/UCjyB1FRMHeuVlMPZFFRuObOoTMUmpOhNQS1q4gMSM8xQLm5uZ6BzwM9TvrP5LocLwpeRPPnz+e5557zdTXEDzW3t3FwUQQAlXcv47rbl/i4RiIi4tafv9/qKRK5QDp7jCPq0JgiEZGAo1DkD8rKwGQythKYysqIjowlxu5k/F+dRO7ar3YVEQkwCkUiF0GnFoQVEQk4CkUiF4FToUhEJOAoFIlcBApFIiKBR6FI5CLQ5TMRkcCjUOQPJk+GPXuMrQSmyZM5+PLfaRwJFXeZqE2OUruKiAQYrX3mD8LDYdw4X9dCzkd4OG3Dk3GGBuEMhU5Tp9pVRCTAqKfIH1RUwH33GVsJTBUVJC/5CeHHnAx9y0n48Rq1q4hIgFEo8gc1NfD008ZWAlNNDbEvvs4m573E2iGkvlntKiISYBSKRC6gevNQY8fZ4duKiIjIOVMoErkIXApFIiIBR6FI5GLoVCgSEQk0CkX+YNgw+P73ja0EpmHD2HXXVVQmjqFkzj00R5rUriIiAUa35PuDYcPg4Yd9XQs5H8OGsf3D0zi6zcrRZCtXhn9P7SoiEmDUU+QP6urg5ZeNrQSmujqSN+whtv4k4/eVEdFQr3YVEQkwCkX+YO9e+NCHjK0Epr17ufYPr2E9vIUH//N9LCcb1K4iIgFGoUjkIgjtbPN1FURE5BwpFIlcBGEdCkUiIoFGoUjkIgjtaPd1FURE5BwpFPmDsDBISzO2EpjCwmiyhNAWEsHJuGGE0IYrJETtKiISQHRLvj+YMkWDcQPdlCm8+5kRbGcu29Pm8kDNXXSmjSZY7SoiEjDUUyRygYS0N3v2XQTR0uTwXWVEROScKRT5g82bITnZ2Epg2ryZqx4/ztTdb/O9wk8QUm2iff17alcRkQCiUOQPOjrg5EljKwGpo7WVkBYXwc52opvrwBVEW32N2lVEJIAoFIlcAFUfFHl97zKZqa074aPaiIjIQGigtcgF0PLSz7y+d2KituaIj2ojIiIDoZ4ikQvAmTDa6/vG4AgaHQpFIiKBRKHIH4wfD++8Y2wlINVPy+CDu+LYO2oGv1vwU6pTEqkJbVK7iogEEF0+8wfR0XDVVb6uhZyHNlqpGB1Pgymehqh4Roel4GqrUbuKiAQQ9RT5g0OH4BvfMLYSkMwnqpn2Vg0jj+3mo288hrk9goTKE2pXEZEAolDkD44fh1/+0thKYKquY3JZLUmOw1y/4VnMnRHE19SqXUVEAohCkcgF0NbW6PW9MzyWqLYmH9VGREQGQqFI5AIwdbR5F0TEEerSpI0iIoFEoUjkAgg+JRS1Bkf5qCYiIjJQCkX+ICkJvvAFYysBKdjcQvn1E6iOHco7V3yY+jGT6QwD+7Xj1K4iIgFCt+T7g9Gj4fe/93UtZIBcTiehoS2U3f8hDmyaxIHhk1g0dw4VaxPgpiuxjh599hcRERGfC7ieIofDQX5+PtnZ2Zf0uRdVUxOUlRlbCTjO5lrMHU4sx1qIbKplxPG9BDU1Y3IGE1lxTO0qIhIgAioU2Ww2li1bxvLly6murr5kz73odu6EzExjKwGns/4EobWQvbSQ8QfK+Orfv07w3j1E1Lq46uevql1FRAJEQIWirKwsCgoKsFgsl/S5ImeydufaXmUulwuX2eyD2oiIyEAFVChyS0hI8MlzRfqy7+C2XmVOpwtnkIbsiYgEkoAMRSL+JKW9iaagEK8yp9OJS6FIRCSg6Lf2OXI6ndTV1XmVhYWFERYWNvAXDQqCmBhjKwHHXHechvAonNFhOE1mWkIj6HCacAUH0x4aRIjaVUQkIOi39Tnas2cPcXFxXl/Lli07vxedMQPq6oytBJzg+hMcGZXMsZ2H2Dzher73hdXUWyfSMNzCy9+YqnYVEQkQ6ik6R+np6axbt86r7Lx6iSTghdcf51hkPB0dnZ6y9nYn5pAw2lvqzvBMERHxJ+opOkdBQUHExsZ6fZ13KNq+HaZMMbYScKLrjhHUFMnQrKuYVP4+3/jbFwjZs5O4mg5u+rNd7SoiEiAUivxBS4vxh7Olxdc1kXPkbGkgsqWO9vB4QvfuIqy9iaHVB3E2NRNKCLHV7biam31dTRER6YfLMhTZ7XZfV0EGifpjuwFISLF6lXd0OAkJCQeg8ZTFYkVExD8FZCg604zUy5cvJy0tjeLi4nN+rsi52tQ1ceOMcXO8yjvanYTEGAvB1jbXX/J6iYjIuQuoUGSz2cjLy8PhcFBWVsby5cspKyvzOsZqtWKxWLBaref8XJFz1Vy5l5agEEakjPUq7+hwEjr+GgCKX/8LT+183we1ExGRc2FyuVwuX1cikMyfP5/nnnvuwr6owwFvvgnXXw9ahiSgvPTzOwip+IB5P97BrsJintltYnrQSVpmX811Hx1D5f2TWTMljZ9Ou41D9z/q6+qKiAxa/fn7rVvy/YHFAvPn+7oWMgANh7fiih0KFguVc26mqWY3h8dMISI0hOTho3n66kyuOVhKiLMDl8uFyWTydZVFROQ0Aury2WXr2DFYtszYSsBwupyMrD8OwybAsWMM/9vvGVGznwzbk5hPHodjx5h7IImE+iZuPLGLls52X1dZRETOQKHIHxw5At/+trGVgOGoPkxCexPJqbPgyBHGPfELhtfsJ/OZPxB6shKOHGHun59hd/Aw5h/ZSG2bplwQEfFnCkUiA1RzYBMA4SOmeMqCzWbAGGjt9nrSBKbUHWHz+8Ucaqjhx+teQEP5RET8j0KRyAA1HdxChymI2BGTPGXBwcY/qZ6hKH/B99kdPYS6fz/M3NWP8qetb+Jobbrk9RURkTNTKBIZoNqK9RyMSCAxOsFTFhzqDkXd66AlRMTw9oyPcUXdYa6qKgegqUPji0RE/I1CkT+wWCAnR7fjBxjHnnfZGTuUlIgYsFgon3kTnSnDOHbNrTQGR3q161fv+Qnr4sfy5b1rCO9so75d44tERPyNQpE/sFqhqMjYSkCorq/C2niCXTFDjdvsrVb+e+/DNGRexfbv/5YTMUO82jU5Ippfp2cR397EAxVvU9/Wwk/LXuGnZa/4+kcREZEuCkX+oK0NDh0ythIQdm19lWCXk8Uf/opR0NZG8LEjRLraiHIcp7O5xatdg4PMHImI5/HU68g5XErh6u/z602v8utNr/r2BxEREQ+FIn+wdSuMGmVsJSD8fc1K2kxmJk6+CQDXli186ef3MGrDG8xZcC2JR+10btrs1a6H7n+U4hGZvJtgJW/Tv0hpqfO8XnNHGyOfeIi1R/b65OcRERGFIpEBmVh3hPLoFMLDIgGor28FICY2zHNMa2tH7yeaTBRMvI1mcwjLtvyLqI4Wtlcf5URzAwBP7frg4ldeRET6pFAkco5cLhczag8RPeE6T1lVlXGLfXx8hKestaWPUATUhUTy0LQcEtsa+OG2Z/nIv3/Om0f2GK+N5i8SEfEVrX0mco5aju4iubWePda5nrKONuMW/JDQ7n9SLX30FL2bs4QOlxNHazPfbm/mp5tX84Ntz/JdgKBg3j1mp9PpxByk/6+IiFxq+s0rco42vPMPOkxBkDbHU9bWFYrckzdC3z1Fo2ISSI1NYmbyKLZYRvKdqR9jpuMAP9n6DOGdbVS1NPKP3etYtWc9Py17RTNfi4hcQuop8gczZkBLC4SE+Lom0g+71z1DbMxQgkO6L5U5xk7gu1/7Dz+65yO0zL+TI99/FcfYiWds169Ov5ngoCCWmoL4ydZn+OXGf/L/pn6MX25aQ2WTMQj7I2OnMjF+KEEm/f9FRORi029afxAUBGFhxlb8WmdnB1OqKthgGUP26Mme8rZ2F+bIcAgOJiw2CoKCaG7pPGO7fivjFr4+I4uffGo5X51xDwltjfx+w1PEVO7xHHOksZbRf/k2z+/bctF/NhGRwU5/hf3B7t1w443GVvxWh7OT7N8/SFxHM9bZdxIR3N0DZK7Yw/1/XwIvv4zpppsY2VSJaU//2nV60kiiUjP5QsYnqQmJ5Hcbnmb+kQ3gcvEfu7Ho7KuHdl3MH01ERFAo8g8NDfDGG8ZW/FZ1SxPXntxDXXA49936Ja/HXLUNjNm3GU6ehDfeIMbcToejrt/tGmQyURUWzVdm3ssLw67ga3tsPLz9WWy73gNg1Z71FO0pvSg/l4iIGBSKRPrJ0dbENSf38m5iGiGh4V6PtbZ1en0fHh5Mc3P/F301d40Zag8K5jfpWXx/8h3MdBzgifV/5roTRi/R198u4mSzgrOIyMWiUCTST7WHtpPadJK1SeN6PdZySgCKjAyjqbH/y7aMjkkA4E83fQKAt5LH89lZ97MjZhg/2P4cP9j2H5Ja65m9epnX814/vBt77clz/VFERKQPCkUi/VT17tM0B4WwPn5sr8eam7wDUGRUKI3nEIqWXXUn/7j1s3x4zBQeveouAE6GxfC9KXfy8OT5TKk9zF8/eJyFFW9h21vKt9b+C9vBHdz3yp+5/pmf8d4x+3n9bCIiolDkH0aPhpUrja34JZfLRcym53kzeTxPf/TLvR6vDEtg8xd/ALNmwcqVBI0dw7GQ+H63a1RIGNcNTyfIFMR9E+dgu/NrxgMmE28mT+DTV36W54bP4JP73yVi+c2cfPMJ7i/5i+f5OS+uoLnDCGGN7a0X4kcWERl0LkkoqqurO/tBg1lSEjz4oLEVv2Tf8grJDccZdv0DzBma2uvx485wau7+BEyYAA8+SPjIoZxwhuN84LMDateJ8UPZ88kfcuWQsQCMTUllw9Wf4f7ZD7A9djhLd73I4+uf4ObjOwhyOQHYXn2M/fVVTHjq+6w5uNPzWkcaHJ55j3pq7mjjjv/9gf31VedcPxGRy9EFD0V1dXX87Gc/49ZbbyUxMRGz2Ux8fDxmsxmz2cytt97Kt7/9bfbt23eh3zpwnTwJjz1mbMXvHKiv5l///A4nQ6NJn31Xr8dbWtoJdlST+uozsGsXPPYYlo4GwhoctP+xcMDtGhEcyp/nfYon5n2Kl+/4Ks98+HO8ufh3xOU9yRdnfoLK8Fj+347/8ed1T5BVuY0DtcfZWX0MgBf3b6W+rQWXy8WVRY+SueoRTjTXe73+9upjlJ44wBPb3xlQ/URELjcXLBTV1dXx+c9/nvj4eJYsWUJJSQk1NTW4XC6vr5KSEh599FHS0tL40Ic+xP79+y9UFQLXgQOwaJGxFb/z0eLlZFVu440R05mWPKrX47W1LcTXn2D0j5fA+vWwaBGWuhPE158g7EufP692tYRFek0SGR4cwoL0WeyIHc7SaTl8fuZ9HIqM59s7XyDpZ9m8/NcvE9XRwj/3rGfS0w+T99rTnuf+fIMNgHZnJ62dHZhNJgA6tZSIiAhwgZb5eOaZZ3jwwQcBuPvuu8nOzmbWrFlYLBYSEhI8x1VXV+NwOLDb7bzyyisUFRVhtVpZsWIFn/3sZy9EVUQuqGNNdVx7eAMRne1k3vHdPpfbqHO09CqLig69FNVjV+ww/t/Uj2FtOE7OoVIeqHibT+97h5eHTuU/w2fywv6tnmOf2vU+ueMyuOP5PxISZOaZD38OAKdCkYgIcAFC0U9/+lMKCwtZuXIld9999xmPjYuLA2DmzJncfffdFBYWYrPZ+NznPkd5eTmPPPLI+VZH5ILZV1fFtcXLeeJIGW8npWNN7N1LBEZP0amiosMuat3ez32Ilw9s43vv/5c7rNN51r6Jn036MCM/vpz1//kx849s5M4jG9gaO5wXhl7B6ykTaDGHcsfzfwSM3qLb//d7ANqcvReuFREZjM4rFD300EOYTCb27t074NfIyspi7969PPTQQ/zsZz/j//7v/86nSiIXTMnB7cyttjOmqZpfpt/C90P6Djp1tS2Eh3v/UwoODiIy8uL1Fo2ItrAwfRbvV+7jx3Pm86VpNzE2NoGI4FC2dXRyz+ZXubqqnNuObub/dr/El8rX8FryJF4aOoXtsSNwdV06A/jH7nXcMmoS6ZYhjI1N7HcdnF0DvLVYrYhcLgYcitasWQPAsmXLznJk/zz66KM89NBDbNy4kRkzZlyQ1wwY0dFwww3GVvzGD97/H7/b/y5bYkcwLuN2z51gp6qpaSYsKd5ow6Sk7raMiaZq6mwSL1K7RoWEUdg12WN8eJSn/LuzP0xUcCi/2GjjjeQJDGmp5dZjW7nt2FY+cmwzlWGxvJY8gddSJrEnOgVMJu5f8zcAfn/DPYSazYy3DOF77z1HXFgEf7jxXsC4i21z1WE+NGYKANf962eAibU537ooP5+IyKU24FA0a9Ys5s2bdyHrwqOPPkptbe0Ffc2AMH48vP66r2shXf6y411ONNeTWbOfyfVHWTIth//c8PHTHl9d1UT4xAlQ8LpRcOutANQPPcrbP/gzd4wffwlq7c3dEfS92R/hh+ue529jr+HJMVczrfYQN53YyYcqt/LxQ+s4GBHP68kTeTtpHHuih/DFN/4BQEpEDMe77lbLHjWZ64aP495XHmdv7QkO3f8oAPvrqy/5zyUicjENOBS5xwf15Xx6e870upctpxPa2yEkBIJ0KcLX/t97z4LLxW/2r2VnzNA+Z7DuqepkE9OmpUBrK5jN0NkJISEEB4GzucVo30vcriaMVDQqJt5T5jKZ2GwZxWbLKH47bh4ZNfu56cRO7jxSxicPvMuJ0GjeSRrH2sRxbLKMgiDj18OX3/wnc4aM9YQkl8uFqcflt7eP7OXa4d1Ln7R0tGOvO8HkhOGX4kcVEblgLspv6kWLFl2Ml718bdwI4eHGVvzCDSd2MbXuCI+PvY6SO79+2uM6O504HM2MqtlvtOGqVZ62HFpZzl33zfZJuw6PMv5zkRQezRWJIzzl37/yIwA4TUGsT0jlpxNu42NXfZFvXrGAN5PHM6fKzvItxTzzzu/5/rZn+fDRzaS01PF+5T7q2owB5XmvPe0ZpA3w8Zcf83rvb7/7H2559je4etzV9vW3VrOt6ojXcU3tbZ5xSSIi/uCC3JIP8LnPfY7S0lKysrJwOBxnPHZQjhuSgBHa2c7n7G+wNnEcpQljmZQw9LTHOmqacTpdWOIiej1mNvuu1y93XCYjo+OZPWQsf7zpXr773nN8e9ZtpMYm8YMPnvc61hIZy5XT7mNvfRW/378Na+MJrq4qZ25VOV/f/QpmXByISGB9wljWx4/l1c42Wszeg8g7nJ3sr6+meG8Zq/eWAvDC/q18ZOw02jo7KNpbxm7HcZ6//UuUnTjAsMg4Zq9exhen3cjSWR+6ZJ+LiMiZXLBQ9Kc//Ync3FyWL18OgNlsJjs7m+zsbO6++27Gjh3rObaqqopnnnmGj33sYxfq7UUuiIb2Vu498D4JbQ38ybqAP3YNMj6dymMNAMQn9hGKgk29yi4Vk8nE1cPSABgTk8jfsu/3PPbTa+7mW2v/5fn+vdyHiAgOoaa1iZcPbMcenYI9OoWnxlxFdHsLMx37mVWzj6tP7uVjh8toNwWxK2Yom+KMS3HbYodz+//+wJaqw151yHvtab4z6zYWpGcCxjQA5bUnmP+/PxDTdSef7eAOhSIR8RsXLBQBFBUVkZmZyYYNG7BYLLzyyiu88sorLFmyBIvF4glJ8+bNw2az6RZ88StOl5OPFn6BPxx8n6dHz+W1xb8lMuTMt9UfO1ZPeHgw0X3MS+TLnqIz+Xj6LG4dPZkP//e3HGpwEBpkBsASGsE3Z2axcNwsriwyBlM3hkTwVvIE3kqeAC4XI5odzKqp4IraQ9x2bAufOPg+nZjYG53C5riRbLKMYmvcSOpCjJD4k/Uv8pP1LwKwvfooNzzzcwDquxatNQcFcbjBQeHWNymvO0lqbCI/mjOfH3zwPz47+RpGxSScWn0RkYvmgoYiMILRggULWL9+PRUVFdhsNl555RXWrFnD6tWrKSoq8hxrsVgUisRv/H37O3xr5wvsj0yg7eYvnjUQAVQeq2fosBivgcdu/hqKTCYTCeFRmLvmFzJ3DQI3mUx8fUYWAC/e/mVu++9vceHit9d/nKaONvLfeYbDkfHEjZzMlo52flR7ghHNDqbXHuSK2oNcd3IPuYeNS2cHI+LZGTOMHbHD2BEzjPLoFDq6wldP26uPMqcrgAG8cRi+OSOLx7avZU/tCZ7Kvp+yEwfJTBl9sT8WEZELH4qsVisFBQUApKamsmjRIs/A6w0bNniFpFmzZl3otw9MU6fCwYOQkuLrmgwqNa1N/G3Hu1w7PJ07nv8DD1S8xT0NJ/hixn08mdm/SzrHjtYzeowFpk402tBi8bRl45ZO/vqDZ/j01KkX9ecYqNTYpNPeVj8tqXtw9l1pMwCYlDAUS2gk1rgk9tdXcU3xTzkcGc/hyHheGHYFAENaaplae5hJ9UeZVHeUG07sItTVSZvJzJ7oFHZ2haTdMUM5HBHvNYmkW3NHe9e2jYUvreSdY3a+MSOLB6dcS2xoOPVtLcSEhgOQ++IKjjQ6WJuz5EJ+NCIySF3wUAScdv6imTNnMnPmTL71LU325iU0FEaO9HUtBpXG9lam/f2HADy16wNmV1dw34H3WJl6HQ/duYQR0ZazvkZnp5PjlQ1cOWeUdxt2TdZoCg+jNjrJeMwP/eHGe9lRffSMx4T26N3JSO7urRkTk8iP595hTF8AfG3GPEzALzeuoTI8jjVDjEVsQ5wdpDUcZ1LdUSbVH2VOtZ27D5cB0BwUQnl0Mnujh7A3OoU90Snsi0riZIsxTuuDyn2e9/vFRhtrj+5lbGwSq/as54l5nyJ79GTePWYH4O+7P+Dh9//HI1fdyV3WGZiDguh0OjGZNOO2iPTfgEPRxbqDbFDemWa3Q34+FBSA1err2gwKT+5837PfXnOIb+94nvcTUvncV1aTFt+/HrsTxxvp7HQydHhMdxt+8Yvw+99DQQGxVYe54i8FcOef/bJdY0PDmTM09bSP//r6hUxJGHbaxz81cQ53pF7BSwe2c5d1Bj8te6XXMe1BweyMHc7O2OH82/2+7c2Ma6gkveE4aQ3HmeE4wO1HNmLGRYcpiP1lT/JQdArlUUZI2heZyImwGN6v3Mf7XUHJPQO325K1zwDwtbdW43Q5WZBujIsaGhnLj+feweaTh/j0pKu8nmM7uANrbDLWuKR+fFoiMhgMOBSlpqaycOFCVq1adcEqc+utt3qNORo0HA4oLoalS31dk0EjuGscTURHGz/Z+gwt5mCWTfww2/sZiAAOHqjBZIIRIy2w/YDRhnfe6WnLWGcL47e8YbRvALo7beYZHw8yBREfHsU942cDcP/kq9npqKSi9iQHGqpZmJ7Jqj2lnuPvHX8lf9/9AXUhEZTFj6Wsx6SYYZ3tpDaeYFzDccY1HCe94TjXn9hNeNditY3mUPZHJrIvKslrezwspnv67i776qq47l8/o7KpjsqmOs+cSqeGos/Y/kpokBn7p38CwOEGB3OKHqX4tsXMHep/IVZELr7zmtF60aJFzJ49m+LiYsaMGTPgSuzbt4+8vDyWLFlCbGzsgF9HpL8igkMJcjn57o7/MrzZwVdn3ENdSOQ5vcaB/Q6GDI3ptRisW0KC8XrNze30vmH/8jMyOp6nb3kAR2sTU//+Qz48Zhr76qq4K20m902YAxiXuUZFx5M39XrPpTf3ZbidscP5zt3f4ZMlTwBgcrkY0lLL2KYqxjSeZGxTFdbGE9x0fCcRTmPcUZM5hP2RSeyPTOBwRDyHIhN4/p1KjkbEwylzKY184iEA7p90NV+dfjMAbc5Oalub+dgLf+JgQw0AtoM7PaHI6XLS1NFO9GkWAxaRy8t5jSnKyjLuVJk5cyaf+9znWLx4sdd8RGezb98+Hn30UYqKilizZs3gu2wml8xX31zFv8o3sOyqO/nkxLm4nJ18fXcJs6sr+PWcB7CHJ3CXdcY5vebBAw5Gjbac9vH4BCMK1VQ3DYpQ5GYJi/SsjzZv1ESvx16/6xvEhIYzJDKWG0eMp7WzgwnxQ/jkhDm8fXQv1w9P57fXf5wvv/lPXCYTxyIsHIuw8F5iGounXEvBtrcxuVyktNYxtisojWk8yeimaq6qKieuo8XzXidDozkYYQwEPxwRb+xHxPPUtjd5Ysc7nuOm/P0HXnVs6mhjXeU+0i0p/HHLm/x+y+s8evVdXDXUSnJEDBtPHGR52StEh4Tx91sf8IxZWnNwJ5kpo2nuaGdYVPdyRS0d7QSZTISajV+326uP8N33nmPVhxbxo3Uv8Pj2tZ7PS0R867wHWmdlZbF+/XoWLFhAQUEBFouFrKwsrFYriYmJWCwWEhISqK6uxuFwUFVVRVlZGevXr8fhcJCVlYXdbh+ca57JJfOv8g0ALH33P7x1eDfjbb/io8c2s3zCbTz+uUL+1MddUGfS1trBsaP1zL369D2klvjuUKRVwAzjLN2XJ8fGJnr2zUFB3DDCWDj3rrQZnjve3L07AHenZbBi29u4TCYqw+OwDB3PKkel1+vHtjczormGkc01jGyqYURzDRPqj3Hz8R1Edhq9S06gKjSao+FxHIuIM7bhlq5tHE/teIe/7XzP63UfeuffhASZGRYZx4GG7jv2/rLjXT40egpbqg7z2Vef9JQ/fcsDpMelMDzaQvqT32O8JYU1dxnLxfx8g433K/dR1dLI49vXnsenKSIX2llDUX8GPlutVtavX8+KFSsoKCjwjAvqa+4W93pIGRkZLF26lLvvvnsA1b7MDB8OjzxibOWiMrlcRiA6uomfjf8Qrwyd2ud5ejYVFdU4nS5S07omF3S34bRpnrYMAV7PeoAOUwxTLuyPMWh8Z9ZtFO8t49uzbmNKYve/D/fdZ998u8hr3FJdSASv3fcjMlc94v1CLhfx7Y2M6gpKQ1vqGNriYHizg4ya/SS1NXoO7TAFcTwsxissVYbHciIshuONMYSERdPetVjukzvf53vv/7dXvT/xyp9Jjojm2Y98ARcudjkqGfXEUixhEZ67+No6O3o975+713HbmKlsrTrMyJh4QkxmSg7u4DvvPcv+Tz/imVPqcvb+sQrS4pJJioj2dVVkEDprKHrkkUe45557uOuuu876YosXL2bx4sWe+YjWrVuHw+GgurqahIQELBYLs2fPJisri5kzzzyIc1AZOlSDrC+it4/sBYzbw/9v10vMO76Dn43/EC8Nmzbg19y7+yQxsWGkpHT94u7Zhldc0X1czmLMbZf/H7KL5fPTbuDz027wfP/23d8i1BzsWfD2O7M+7BWKAM+ludcP72bDx7+D7eBOvrX2X9SERvP5q3N5pGuG7Z5CO9sZ0lrHsOZahrbUMqzF2I5rqOTak7u9LssBVIdEcrwrKLm/jnd9nQyL4WRoNCeaG/jca097nuPCRU1rEx1di+D+d98Wz2O/2GDjFxttAJQc3MHLB7YTExLmmfkboK69hfiwSFwuF2P/+h1+ff0C7uy65PubTa8yd6iVK4eM5URzPUnh0X2G/Q5nJ22dnV4Tk9a1tRASZCYiOOTMjXGJ3P1iIRMsQzw9ayKX0llD0aOPPsrs2bP7FYrc3PMRST85HPDmm3D99cbkf3JBffzlx4hpb+YH255lct0RfjTpdvalXQP1VQN+zb17qkhPT+r+w+Nuw+nTYdMmoy2BKw6sY61LszFfKD0vuQEkhEfxuxs+zo0jxrOt6gitzk4AvjkzG0drM/FhkdwzfjaxoeHkvfY0t4+dxvXDxxFkCuLJne/x5K73SY4wAszByEQORhqv/1DmraQlDOdTXYO+wzvbSG6tJ6WlnuRW4yultZ7k1jpm1hwgubWO6M42T72cQHVoNCfCYqgKjaI6NIqq0Giqw6KoqionPTSKwrVFBIVG4jQFeQIRwMsHtgN4BSKAmpYmSg5sJzIkjE6Xky+98U+K95bxZPb9LO+aDsF259fI+s+vWJr5IZo62nhwyrW8emgXNwxPJyE8ksxVj1DV0ugZw/TCvq0sfu0pr7vw3A7WV7O1+gi3jTm/yUe3Vx/llmd/zXu5+YyMju/XcyrqTp7Xe4oMlMnlvp51BkFBQdhsNm6++eZLUSe/Nn/+fJ577rkL+6JlZZCZCaWlkJFxYV97kPr2u/+hrq2F1NhEnnv77/xg27NEdbTx/6beyXXX3Mt3Zt3G99//L/a6kzx1ywPn9NpNjW384LuvkLNwOrPnjDIK3W341FNw331GWwJkZvLre3/JopVfIDLSPydxHCxqWhqJD4/qVVbd2sSjpS8xKX4Yv9hoY0hEDKUf/w4Aj29fy/ff/y/Lr/kY4+JS+NgLfzrt60d2tJLSWk9Sj8CU3NpAYlsDCW2NJLY2YGlvome/oRNwhERSFRbdHZxCo7qClLFfExqJIySSJnNor+kHziY6JIyGU8IV4AlFo55YigvjT0DhTZ/gI2O7e0+vLlrOgYbqXoPAXzu0i+FRFlJjEwkymSivPcm/yo1LnH351cY1/GxDCQBlC79DSmTMGes88omHCDYFse8zj/T5+InmetYeLff0kp3JV95cRWiQmZ9dm3PWY+Xy15+/3/0aaG2xWLDb7X2GosTERGbNmkVGRgYLFy7UHWTiU0caHCRFRBsDZV0uPnxsC7/bu4YDkQl864oFHI2wMLuzHZPJxA/nzh/Qe+zdcxKXC8aN7/+kf4cP1ZF+DsfLhXdqIHKXxYdHsfLmTwLwjZlZXo9/dvI13Dv+SiKCQ3B2XfbqS2iQmabgMPYFh7Evyrudo4JDaewwepGCXE7i25pIaGsgsa3RCEytjZ7glNp4ksyafSS0NRJyyvu1mcw4ugKSIySix36k135NaCS1IRG0mEP7DETgPYDd7Zcb19Dc0UZcaAS2gzs9A8o7nU4ONdaw8cQhvvjGP077GXxzZjZh5u4/KduqjvBcxWbPnGAARxodZwxFa7sudXe4nKzes54F6d5LQXU4O8n+z6852dLADSPGc7K5gZ+sf4FZKWN5cPI1hJ9yCfCZrhssBhKKjjfVnzXAnS/3rOufsf2VW0ZN5r6Jcy7q+8nZ9SsUORwOHKeZgK6mxpjbY8GCBbqDTC65DmcnJkyYg4JoaG/lyqJHuW74OBJaG/jm7le4qrqc54dewW/S53kGx445z5XXt245xrDhMcTH9+9G++BgM4cP1SoUBSj3WJsgUxC/vDaXr79dxK+uW8AVSSPYdOIQz+/fwicnzGVsbCI3PPNzABLDo6hqMQZvD4mMxd51OchpCuLOmbfy+Pa17DnDe5pcLmI6mkloa8TS1oylvYn49ibi2po8+0Oba5lUdxRLexMxHb3DT3NQCI6ugFQTEkldSAR1weHUh0RQFxJObXAE9SHh1AVHUBcSwb6TbXztzdW9eqPG/PXb/fqcqloaeXH/VoZFxvHhsVP5v7X/YkvVYa9jnt79AamxSTyw5q/8/sZ7GRrpPS/dwpcf8+z/aeubvULRrza96lkGxr1MDxhzS71ztJynbrm/z2Vdbv73L5g3chLfmX0bLR3tjHvyu/w16zO9poxwW3NwJ5+2/YVX7vgKkxMu3g0wE5/+vmetv1cP7eK+iXM41lSH7cAOnLhYmD7LK2jKxXfWT7u2tvaMj5tMJoqKijTpovjE7NXLONHcgDU2iXZnJ7hchGx4jj/vXUOHKYj/N+Uu3kka5zneGpvEZydfM+D36+hwsnP7ca67of8zHienRHP40Jn/HUlgyE3PJDc90/P9eMsQr+//dVse8eGRjLcM4WRzA62dHXxQuY8vv/lPbhgxnjcO72bBuExSImJYVvpSr9d/9iNf4I7n/4DLZCJ74jXGVBJREBcaTmbKGP59aFef9Qp2dhLX3oyl/fQhakRzDRM7WohtbyG2vRkzvUdOtJnM1IWEUxcSQX1weFeQMkJUXdf37iDVEBxGQ3A4DcFhNJlD+efudZ6xUT+eewe7T5kuAeAfu9dhjU3i/cp9zFr1CEsybuEr02+m5MB2piaO8Dq29ZS787ZVHWHDiYOnbZs3j+xhxda3+dy066ltbSYurPs/Lbsdx9ntOM5nJl3Fvq6xhE/ueo/goCCu6xpj1tzRzm5HJdOTRnrqvttxnFue/Q0/vzaHhemzcLqcLC97hfmpVzAqOgGzKYjIkFB21hwjPiySiOBQ/luxmfx3nuGZD3+OK4eM9arj1L//kDut0/nx3DuA7sWPe/rKG//kna41/Y431/O16TcT3GMNwv74oHIfIUFmZiaP6tfxT+58zxhOcMq4ssHorKFo/fr1gHHb/emcGogG5fpl5yM8HCZPNrZyTk40G/9rtNedxNpwnJ/vfZWZtQd5NXkiv0nPoi7E+MV47bBxvH10L5+aOPe8Fggt33uSlpYOpl4x1PsBdxvGxHi35eTJpIxJZG1FNS6Xa0C3/0vg6LmWnPuW8ruiu+ddqmltIj4skimJw7nTOoNDDTWkxSVT19ZCfXsL05NGMsEyhHRLCmlxyYARlCYnDKOurZnMVY8wOWEY27sW8r16qJV3jtnpCDJTFRZNVZj3bexv3f1/XPevn/WqZ2xIGHHODjoaqojpMEJSbHszse79jhZi2puJbW/B2nqiq6y5zx4pgE5MNL7zOzKCw2gMDqN+0yqW9ghNDcFhNHbtv/rSXqZ17f/53X+x5chuXjxa0auHan99NZ94+XEmJQxjw4kDnnXvzuRXm9YQZDLxw3XP990+Rd3jo2wHd2I7uNNYt++au8l5sZBNJw8xPWkks1OM+cdqW5sB+ObbxTS0tzItcQS/2/w6VS2N/GP3OoZFxrFu4VKy/vOrXu/1sRf+xFen38zXZswjpCvUOFqb+MuOd/nc1Ov7DEQATT3Kf7VxDRW1J/n9jfd4HTPyiYfIm3Id373yI17lzR1ttDudnrFv/Z0U9NebXqWt6yaFgahpbcKEMXFroDtrKCopKcFkMpFxmgHAfV0yW7duHbm5ueTm5rJgwQIFpLOZPBm2bfN1Lfzausp9BJlMZHb9stpadZiorqUXElob+NT+d/jI0c0ciownf1oO6xJSyUwezR9uvJekiGhCgoKoa2shLvT85pbevPEoCYmRDB12yliDnm04v8dYpW3bSNlaSe3j66ipbiYhMfB/acjAxff4ozEi2sKIaAuA15w87lvR252dWOOSyUwx7l6MCA7x/JE70uCg3dXJmJhEalqbqGyqY7wlhev/9XP21Vfxwu1f4khjLamxSUyKH8qOmmNe9XhwynWUnjjAG4c7OdpV9otrc6hqaeQnfUxZ4BbkchLd1dsU3dFKdEcLUR2tnv3ojlaiO1sZHmQirNFBYutJojvdj7cS5uw9NxPvwTcw0RAcTmNwKE3mUJqCw2g0h9K0w9jOMIeSHhxKkzmMpuBQGs2hNJtDaezqpfJsXa7TBqLT+bd9I1EhYWw6eQiATScPefZ3O457jvv++//l8a6xZ01dY8SONtUy858/Pu1r/3rTqzxTvoG8qdczrcc8W3OLCk77nFP/3/RsxSZ+f+M9HGuq4+0je8kZZ/wtXrn9bb575Uf47abXGBkdz11pM5j69x/26mFzO/U/Za8e2sWk+KEMi4ojqKu80+nEHBTEoYYahkTGEhJkpq2zg1BzMK6uz/aTE+Zg7Qrsbu5LmZfDzOxnDUVlZWVYLJbTLt+RkNB7fMaiRYuIj49nwYIFFBYWkpWVxcKFC/nYxz523hWWwemuHv/zOdncwIee+y1JrfV8+cD7fOToZlrNwfwh7SaeGz6Dzq7/lc0aMtbzRwfO/38xba0dbNp4hOtvtJ5Tj89YazwmE9jLqxSKpN9Cgsx8dGzfc2kN73Fex4dFesKWNS6JffVVpFtSuCJpJAAvzf8Kaw7t5IE1f/M8J9QcjCXM+z8IueMyMZlMrKvcxysHd3jKJ1iGsKvrcpLTFERdSCR1IZFMsAxhR0MNjR1tTIwfys6u4LXzvh8Qbg5m7F+/0/tncnZ4Byiv/RaiOtqI7GwlsrONyI42EtoaGdlRbXzfVeZe9+50mswhNJnDugJWV7jqEbaag0JoNofSbA6hxWzs76kq5wpzaNf3RllLUAhPbX8bely6OtToAOBZ+yZPmbu3+nQONtR41vk7m6l//yGO1qZe5c0d7fxo3fM8a9/kCT1Ol8trwHzx3tJegehgfTWxoRG8dXQvn3vtaXZ84mHWHNrJ8CgLnyp5gnFxybz+sW96QlFLZzsRphDmFhXw2cnX8PH02WQ/+yuKb1vM1MQRrNz2NhtPHOTRqz/Go6Uv8fCcj/LCvq39+tkCxVlDkd1u96xx1pfT/XHIyTFG+1dUVGi80dls3GjMa/Pmm6BetTO65dlf03BwK187XMqHjm2l2RzCk2Ou4j8jMmgM7l608/GbP8n1XctGXChbNh+jrbWTWbP7uE7vbsOVK2HRIqMtAa6/nsg332TosFjs5VXMurJ/1/hFBuLX1y3gxQPbiAjunv7BHBTELaMn851Zt3l6gcLMZn489w6uGZbGT9a9QG1bi+d3+Y/n3sGM5FEsL3uFkCAz/7v9i9S1tRiLKJtM7K+rYsOJg547pV7av43xlhTKa0+wperIGRfPXXfvw8z454+pCY3it9d/nK+/tdozmeXwqDicLhfHmurO+DMGuZxEdrQR1dlKRGdbV5BqI7KjlaiurTtEeUJWRxvxbU1EdbYS3tlORGc7EZ1thPfVc3WKNpO5O0Ct+zN/MIfQHBRCS1ewaj4lTPUMW8ZxxldrULDXti0oGOcpl/L7CkQA6U9+17Of/84zfR7zxpHeQ/evKl4OGGPSACY9/bDX43trT/CHLW9wqMEBwIf/+zs+0jUv1ePb13qWoSk7cZBR0UYHyLrj+5n3n18CxkSjl5t+haLc3NzTPl5eXs5rr73GTTfd1Oux+Ph4BaL+cDqhvt7YDnLvHbOT/86/ae1s51CDg6/NmIftwA6CXE6urtrLnZtWkeE4QHVIJH8dew3/GT6T5q4/AIfuf5Ta1mbanZ0XZYmAdR8cJG1cYt+9Pe427OjwbsuufWtaAju2He/9PJELKD48invHX9nnY5+fdgMT4ofyqZInGBoZR3xYpDHdgDmUldve8hw3PNrCV6bfzKyUMaRExBARHOoVsqYkDvdacuVDY4xFbKxxyWSPnuwp3/jx/0dwUBBNHe38acsbfDT1CpIiovn9DfeQFBHNNcPSuGnkeP6+ex2fmXgVJpOJr7+1mv/1mOnbzf6pH2P92/8DjN6qhpBwGkJ6j8Hs2WPVH0EuJ2GdHUR0tvUISt2hKaKznXDPY70fj29vYlhLba9jQ139G5/TZjL3CkotQcG0mkNoDQqhxRzc57bVHOI5riXo9I+1BgXTFhRMbVvLaevQc4b38toT/Gbza72O+euOd0m4DMYL9cdZQ5HFYiEtLe2Mx2RlZWGxWFiwYAHZ2dm6TCb98px9EzeNnEBMaDhFe0rZePIQe2uPU157wjjA5eK/b/+drModfO/4DpLaGtgaO5wfT/wI3/j0L3hkSCrfbm/1ujU3LuzirEd/4kQD9r1VLLx3xoCeP258Emvf2seJEw0kJ2tNJ/GNm0aM52/Z93NTj17Ungvw9nT1sDP/3j8b939MLGF4zQl2h3W6Z98SFskXeizjktA1l1TWqInMShnD9KSRxIZGEGoOZn7qFTxXsfm077fnkz8iyGQirSs89WVh+ixW7Vnv+d5pCqI5OJTm4FDCzSG0dJ750lx/mZ2dPcJTO6HOdsI7OwjrsQ3r7CDc2U5YZzvhzo7ubY/HYtpbSHbWE+bsILyznTCv487ey+XWZjLTajYCkjsoGV/mU74P7vO41qBgXt/zOtmnOdZ9zP6TBxiTFNgz+J81FFmtVsrLy894jMvloqamhhUrVrBixQoAzyW3/fv3M2ZM3yuJ6y61wau+rYUvvPEPwszB/GDO7Tz0zr8BmJMyFmvDca6uKuem4ztJbTpJTUgkr6VM5KUhU9kbMwSAPySNIiTITHxYpOc26Itp7ZsVREWHcsX0YQN6fnp6EsHBQezYdpzkGxWKxDdMJhM3j5zg62qclnts1P2TruaGUy5//+HGe3lwyrW0dXbQ1tnJva887nns2Y98vtfabSOjLZ7LQsYxX+Dvuz/wfF+68Ns0trfR0tlOamwiEcGhvSa1vH54Om/2cVnKbf3Cb/PJV/7cayB7Z5CZxiAzjcHdvVljYhLZdoalhcbFJbPX/R/CU9w2Zgov7ve+Gcfkcp0SroygZASvdk+QCnW6yzsIcXYQ1tnhKQt1dnr2I7t6vkKdXY93HRfadVyYs52z3bdbMzqNMR/91lmO8m9nDUVZWVnY7fbTPm6xWFi6dCnl5eXYbDbPse671qxWq1cvUlZWlueSWmFhIX/84x8v0I8ivvLYtrd595idx+d9qs/H3zlajjUumf9VbKbd2cnnp91AZde4gdbODr77dhGzHAe5uqqcq94rZ0hrHY3mUN5NTKPQegOl8WPoDDIzNWE4t0TF8crBHV7jFnreBn0xNDW2se6DQ9xwk5WQ0HObL8QtNCyYcelJbN9WyfU39n+OI5HB5NbRk/n1plexxvY90WlGstELsb36iKdsetJITzlAuDmET02cw+Kp1zNr1SM9njuKxPAoT09RmDmYIZFnHt6x7Oo7Kdz6FjeOGE9SRDTz//cHAL50xY1sOnmIlIho7k7L4MfrXzjrz9bWYxD0v27L4+4XC70efzL7fixhkV7jfr42Yx6/2riGeSMneoUi97QMLeZQWsyhnGkWtPiwSGq6xipFh4SRMy6DP+5496z17cXlItjl7ApT7UZQ6hGwpsQmUHDNvef+un7mrKEoLy+P7Ozs0z6+YMECvvWt7mRYUVGBzWajqKgIm82YyOvUXqSMjAysVivFxcUKRQATJxprZU3se3ZVX3n4/f+yv76aR666k1CzmbjQCJ4p3+C5SwWg9PgBHv7gf4AxSNASFonT5aS1s5M3D+9m3fH9/Gnrm16v+3x5GZ0VpXzScYAZjgNMqTtCqKuTY2GxrE0ax7sJaWyyjKIjyMztY6/gsbQZ3L/mb9w2ZgpfvOLG087vcbG89+4BXE4XV18z9vQHudtw9GjvtuyxP2nKEJ59ZitNTW1aB02kD1ckjezXbd2TE4bz1C0PcN8rf2bFTfd53fCz91M/AroHLSeFR7MwfRYmk4mxsYnckTqdZys29TlT9KoPLeJQQw1vHdnLhhMHGBOTyCNX3el5/NmPfJ599dXcnda94PmHx07ld1tep+SOr1Ld0sjnXv87n5441/N78UOjp3C40cEDk67mV5te5W/Zn/HMQdVTYng0kSHevxe+MWMeX5x2I+HmYP5v7b885Q9OvoZvrX2GzjMsPeO25d7vsWLrW/xw3fM8evVd3JE6nUMNNWSPmkxlUx1/3vGO57P6f7M+zP2Tr+a2537Dbsdx7k6byQeV+7h+eDqJ4VFEhoRSXnuCor1lvd5nXOp0QuJH9CoPNP1aEHbNmjVUVFTw4IMP9nqsoqKC1NTT/099zZo1lJSUUFxc7NXj5D6JOzsHPmGUL1yUBWF9aOQTD3HfhDk8evVdfT7mFm4O4Qdzbif/nWe4f9LVTE8aQc64zD7XUJqdMoZ1x/cb37hcpLTWM6n+KBPrjjKx/igT6o8R7uygwRzGJssoNlpGUWYZQ0VUEphMni7rifFDeXn+VzAHBVHZVEdyRPR5Tbw4EK2tHTz6o1eZesVQ7l5wxXm9Vq2jmUd+uIbcj0/XXWgiF1mHs5PxT36Px+d9ipt6XDKsa2vhjcO7uT31/P49n437d+PpQt6pvzvdx7nLV936INcM756Nf+3RckKDzLx3rIIHp1xLRHCI59ghkbGULvx2r9d8+pYHel2GPNXiV5/ihf1be9X13+UbyRo1kZjQ3gPa+/q9/5Gx0yi86RNnfC9fu2ALws6bN++0j50pELmfO2/ePB599FFqa2tZvXq1JyRpdt8uBw5AQQHk5xs9DZeIvdZYj+mpXe8zOiaBl/Zvo7G9lT/edC/jLUO8jm3pbPdcj39ixzsA3Dp6itcxoZ3tjGmqxnJsK4sbT5LadJL0+koS2o3/hRwLi2Vn7DCeGHstGy2jKI9OwWkK4pszs7guOp6ccd3LJRxrqmNIRIznHDlbN/fF8vabFbS0dHBzdvqZD3S34X33wVNPGW0JXu0aZ4nAmpZI2frDCkUiF1lwkLnPZStiQ8MveiAC43dW5VmmFwD44ZzbKT1xwKvMhMkrEAFc0zXwffYpS4eAMT7KbUbSKDaePMiHRk85ayACY7HcF/ZvZWSPua+APgff92XVrQ+y8OXH6DiPGbH9ySVdaS4uLo5FixaxaNEiHA4Hs2bNOvuTBoOTJ+EPf4DPfvaChqKtVYcpPX6AjScP8svrFng9Zq89yfXPdE//3/O2zO+//1/WHu09uH7TyUPgchHf3sTwZgdf+UUun2mqYWzTSVIbTzK82eFZT+loeBz7IpN4YdgV7IgZxs7YodSEGgOMPz/1Bv4y8Ur+Xb6Rr06/uc9wfOpCkb7Q3NzOm6/ZmXv16LMv/upuw6uv7m5L6NWuGbNGULxqMw5HMxbLxblTTkR875U7voKja5mQM3lg8jU8QPd6jFvu+W7vaa1Po2yh9wSZH+QaPThXFj1KRz+nBYgNDWf7Jx4m+Bx64cfEJDIrZTS/vn4h1V0LH7crFJ0fi8VCvvt/03JWJ5sbiA0NJ9QcTE1rE9P/8SNeuP1LvRZRLNpTagxq3reZldve9pTvq6ui+LY8NlcdZlxcMt99v+8uxCCXk20VG0lrayCptYGU1jqGtdQyvNnBsBYHw5trvWaUrQ6JpCIqiQ8SrFREJVERlcS+yCTP3EF9+c7s2wBjEKE/K3lpN52dTm7OGnf2g/tp6hXD+Pe/trKx9DA3zrtwrysi/iUxPJrE8HO/0zS+a1qC/kiJ9F5uaHi0xTOg+5702f1+ndg+LpGdydqc7nHEsaHhhJmDeeA8Ftr2Jz4LRWAsByK9HW5w0NTRxtjYRBa9+iQT44fyu82vM2/kRP54473c98qfcbpcFO0to7KpnrS4ZCrqTnK0qZYla/ue7bS0soLpK79KQlsjSa0NJLY1cF9X8ElsaySxaz++rdFr9ewOUxDHwmM5Em5ha9xIXhkylSMRFo6Gx3EkwkKLue/w8+Dka3hwyrU8/P7/WDTlWuLCIvpcNNFfHTtaxztv7+PWD08gJvbCLdQbERHC5ClDWL/uEDfcnKZLyCKDlHsG7wst1Bx8SdcgCw4yU/6p06//Fmh8GorEW0tHO7/fUMIvN64Buq8N2w7uBGDNoZ1MePK7RHe0MLyjhbXvFbGlvcVYvbpr++X2lu4Vrj2rX7cQ3em9urUTcIREGitrh0ZTHpXC+wmpVIUa35/sKneERvaaih6Ma94uXHx64lX8dadxe+e/bsuj3dnJx19+jDBzCCOj43ls3icv7od2EbhcLp59ZhsJiZFcd8OFv91/7tVjWPGH99i7p4r08X3feiwil7d3c3SlxB8FXChyOBwsW7aMsrIySkpK+v08u91Ofn4+VqsVh8OBxWKhoOD0KxVfSnWmJiquGsGzq75OS7SZr7Q3E9W1SOJnOlqJ6mwluqONqI5WIjtb+5xAqzUomLrgcOpCwqkLjqA+JJzj4bHUB4dTFxLR9VgEVaHREDeEPZ0uz8Kpp8odl8G7fdxy+ZO5dzAyOp7fbn6NH8y5nf/YN/L9Kz/KX3e+y8L0WcwZmkqn08lXp9/Mg1Ou7fX8T06Yw7Qk/79ls3TdIcr3VvHZxVcSHNzPeYlSUuDrX4fx441tSopR3nO/S9q4RIYOi2HtWxUKRSKDlDno0t5JK/3Tr1vy/YXNZqOkpITly5eTkZFBaWlpv55nt9vJzMxk5cqVnoVq3eu5FRUVnVMdLsYt+aXvrCJmhTHpVVVoFI6QSBqCw2gIDqcxOJRGcxiNwcZXgzmM+q7gUxdiBJ0hiaPY1mPm1r58YdoN/K9iC1+fMY/c9MzTrsYMxirXZlOQ1yKEGz7+HZIjYvo8vt3ZidlkuuS3y18MjppmfrH8DaZMHcrCT8y4aO/z/nsHeGb1ZpZ8+yYSk/o/hkBERAamP3+/A+qvWFZWFgUFBVgslnN6Xn5+PgkJCZ5ABFBQUEBxcTHFxcUXuJbnLr25kfDjLnbOuIc/f/z3LJr1Gb4+4x6mLXmFRyd+hNT7/8SGK+9l1agribr+Ad5MnsDG+NE8cve3+fa8+/n77V/i/klXA/DtWbdx8DPLPN9fMyyN8k/9mG/Puo13cpeQm27c9n710O5ZlacnjQQgtKvnKMIc0mvK/Og+Fl90CwkyXxaByOl0UbxqM2Fhwcz/2JSzP6GnhgZ49104dszYNjR0lzU09Do8I2MEkZGhvP7qmZfQERGRSyfgLp8BJCQk9PtYh8NBcXExixcv9ip3Lz9SWFjoFZZ8wVSxn9Evuoi/cyy333QfFXVVpMYmYjKZPAPmFqR3T19wp3UGf935Ltf2mMfiR3Pnc6d1OjOSRmEymfjR3PlMTxrJTSPH9zlz66+vX8hDjbdiMpkYEWXhtUO7SImM5c/b13q6de2f+jFNHW2sPVreKyRdjt58rZzdu07w2bwriYg4x593927jdvynnjLmKnL3Yl59tbGfkeF1eEiomRtusvLyi7u4OWsc8QmDYwVqERF/FpCh6Fy4lxpJS+u94rPVavU83l9Op5O6Ou8JucLCwggLCzvNM/qhq5fFFBJurBcXd+ZxJtePSOf6Eb0nE8xM8V54N2dcRq9j3CKCQ7D2mGr+1jFGz8hvb/i4pyzUHEyoOZiPjJ129p8hwJXvreKlF3ZxU9Y4JkxMOfsTLoCrrh3LG6/ZedW297xnyxYRkfMX+Nc8zsK9tEhfl9zcPU4Oh6Pfr7dnzx7i4uK8vpYtW3ZedYyeZszbEzX1lvN6HRmYutoW/v63MsZaE7jlQ2efAfZCCQsL5oabrKx7/yDHK3tfYhMRkUvrsg9Fbn1dcnMHpZ5rsp1Neno6tbW1Xl9Lly49r7qZgo0OO1Mfl7nk4mpt7eCJx9YRFGTiE5+cidl8af9JXHNdKnGWcP737PZL+r4iItLbZR+KrFZjQHF1dXWvx9w9ROcycDsoKIjY2Fivr/O6dAYQHAxJScZWLpnOTidP/62MEycauH/Rlec3SaO7DcPCutuyH+0aEmrmo/Mns3PHcXbtOD7w9xcRkfN22f8Vdoeivi6RuYOS+xifueIKOHHCt3UYZNwTNO7eeYL7F81m+IjzXGutZxv2HLjfj3adesVQrGkJPPvvbXx9XCIhIf2cG0lERC6oQdNTVF7e+9Znu91ORsbpByPL5cnlcvHf/2znvXf2c/eCaZdsYPXpmEwm7sqZRk11MyUv7fZpXUREBrPLMhT1HCNksVjIycnpdZeZ3W7H4XCc93igC2LbNhg3ztjKReVyuXjxfzt5+80K7sqZyuw5oy/MC7vb8Nlnu9vyHNp1yNAYsm9N543Xyjl4wHFh6iQiIuckIENRX+OD3JYvX05aWprXpIwFBQVUV1d7BaP8/HyysrJ8PkcRAK2tUF5ubOWicTqNHqLXXy1n/l1TuOqasRfuxd1t2NDQ3Zbn2K433JzG8BFxrHp6I62tHReubiIi0i8BFYpsNht5eXk4HA7KyspYvnw5ZWXea3S5J2XsOU7IarVSWlpKYWEh+fn55ObmYrVaz2ntNAlsHR1OVv19I2vfMnqIrr3+wi/0er7M5iDuuW8mDkcz/y7eQgCtwCMiclkIqIHWWVlZZGVlUVhYeNpjcnJy+uz9sVqt57zOmVweWlraefqvZezdc5JPfCqDK2YM93WVTitlSDQfy53GP5/eSNq4JGbPGeXrKomIDBoBFYpEzlXVyUb+8tg6HI4WHlh0JekTks/+JB/LmDWS8r1V/PtfWxgyNJrRY+J9XSURkUEhoC6fXbbGjYOXXjK2csHs2XWC3/zybTqdLr709WsubiByt+ENN3S35Xm06513T2XEiDj+8vh6aqqbLkKFRUTkVApF/iA2Fm691djKeXM6XbxasofHV3zAqNEWvvS1axgyJObivqm7DUeO7G7L82jXkBAzn/7sLEJDzDzx2DpaWtovQqVFRKQnhSJ/cPQoPPywsZXzUlvbwso/vsfLL+7ixnlpPLDoSiIjQy/+G7vbcOPG7rY8z3aNjg7j/kWzcdQ085fH1tHW1nkBKywiIqdSKPIHR4/CD36gUHSetm+t5Fc/fZPjxxtY9Pm5fOjDEwkKMl2aN3e34bZt3W15Adp1yNAYHlh8JYcO1vK3P6+no0PBSETkYlEokoDX2NDGP57awF8eX8foMRa+8a0bGJee5OtqXTBjUxP4zIOzsdurePIvZbS3KxiJiFwMCkUSsFwuF5s3HuFnBa+zc8dxFt47g888OJuo6EtwuewSG5eexKfvn8We3Sd4YuUHtLRockcRkQtNoUgC0skTxq32T/21DKs1gf/Lv4HM2SMxmS7R5TIfmDAphc8unsPBA7Ws/OO7NDa2+bpKIiKXFYUifxAfD5/4hLGVM2pt7eDF53fy84I3OHq0nk/en8kn759FTGy4byvmbsNRo7rb8iK0a9q4RPK+eBXV1c384TdrOXGi4YK9tojIYGdyaS2BczJ//nyee+45X1dj0HE6XWwoPcxLz++ksamNG29O48abxxEaavZ11XzixIkG/vLYOhoa2rjvUxkBMSmliIgv9efvt2a09gctLXDokDHHTbiPezz8jMvlYvu2Sl5+YRfHjtYz9YqhfPSOySQkRPq6at7cbZiUBCdPGm0JF61dk5Oj+dLXruXpv5Xx+IoPuP2OyVx93djL+vKhiMjFpstn/mD7dkhPN7biUb63ij/85h3++vh6oqJC+eJXr+FT98/yv0AE3W34/PPdbXmR2zUiIoT7H5zN1deO5dl/b+Ppv5bR3KxJHkVEBko9ReJXXC4Xu3eeYE3JXvZVVDNiZBwPfm4O6eOT1AvSB7M5iPl3TWFsajzFqzbz65+9xb2fmqn10kREBkChSPyC0+li25ZjvGrby+FDtYwabeHTD8xi0pQhl24CxgB2xYzhjBxl4e9PlvGH37xD9ofGc+PNaZjN6gwWEekvhSLxqdbWDsrWH2Ltm/s4fryBcelJLP7CXNLGJapn6BwlJEby+S9fzSsv7eaVF3exdfMxFtwznWHDtaaeiEh/KBSJT1SdbOSdtftZ994B2to6mTJ1CAvuna7LPufJbA7ito9MZOq0oRT9cxO//vlb3Jw9jpuzxhEcPDjv1BMR6S+FIn+QkQGDYGYEp9PF3t0neWftPnZsqyQiIoSrrhnD3GvGEh8f4evqnZ+ebfiJT3SX+6hdR4228JVvXMuakr28WrKXjWVHmH/XFCZOSvFJfUREAoFCkVx0NdVNrPvgIOs/OISjppmhw2K4e8EVzMwYQcggnWfoUggONnPrbROYPmM4z/17K39e8QGTpw7h9jsmk5gU5evqiYj4HY3C9Ae7dsFVVxnby0RHRycbNxxh5R/f49Efv8qbr9sZPyGJL371Gr7+reu5cu7oyysQudvwxRe729JP2nXosBgWfX4u9306g8OHavl5wRs8/9x2mrRMiIiIF/UU+YPGRnjvPWMbwJxOF+V7q9hYdpgtm47S0tLB2NR4chZO54oZwwgLu4xPN3cbVld7t6WftKvJZOKKGcOZOCmF118t583X7bz/7gFunDeOa69PHbQzg4uI9HQZ/5WSS8HlcnHwgIONZUfYtOEI9fWtJCRGcs11qczMHEHKkGhfV1F6CA0L5pbbJnDVtWNZ88oeSl7axdq3KpiXnc7sOaMICVE4EpHBS6FIzpnTaQShbVuOsWXzUapONhETE8b0mcOZkTGcUaMtup3ez8XEhHHn3VO57oZUXnlpN88+s5U1r+zh+puszL16zOXdqycichr6zSf90tHhxL63im1bj7F1yzHq61qJig5lytQhfCz3CtLGJWqSxQCUmBTFPffNJOvWdF5fU86L/9vJa2v2cu31qVx97VgiI0N9XUURkUtGocgfjB0LTz5pbP1IU1Mbe3adZPu2SnZsq6SlpYP4hAhmzBzO1GlDGZOaoCDk5m7DuXO929IP27UvycnR5H58Olm3pPPGa+W8WrKX19eUkzl7JNdcn0pKii6Disjlz+RyDYIJci6g+fPn89xzz/m6GheFy+XiyOE6du04zs4dxzmw34HT6WLosBimThvK1CuGMmx4rC6NDQL1dS28u3Y/772zn4aGNiZMTObaG1IZPyFZ7S8iAak/f7/VU+QPTpyA1athwQJITr6kb93U1Mbe3SfZueM4u3acoL6+ldAwM+njk7jz7qlMmJQS+BMrXgruNrz5Znj1VaMtwWfter5iYsO55bYJ3JQ1jk0bjvD2mxU8XvgBySlRzLlqDBmzRhAdHebraoqIXFAKRf7g4EH40peMOW0u8h/P1tYO9tmr2bvnJHv3VHHkcC0uFwwZGsPMWSOYOCmFsakJBAdrCqtz4m7Dp57qbku4ZO16sYSEmJl15SgyZ49kn72ad9bu58X/7eTF/+1gytShXHnVaMalJ+kyqohcFhSKLnPt7Z0c2F/D3j1VlO856bkkFhsXxrhxSVx97VjGjU9Sb5CckclkIjUtkdS0RBob2ihdf4h17x3gsT+9T3x8BLPmjCIjc4RmyhaRgKZQdJlpaWnnwD4HFRXV7LNXs39/DR3tTiKjQkgbl8QdH5tCWnoSyclRGhsiAxIVHcr1N1q57oZUDux38P67B3jj1XJKXtrN6DEWZmaOYPqM4UTH6PKaiAQWhaIA53A0s89ezb6KavZV1HD0SB0uF0RFhTI2NZ4PfXgi49KTGDosRpc45IIymUyMGRvPmLHx3PmxKWzbVsnG0sP89z/b+e9/tpM+PokZGSOYMm0o4eH6VSMi/k+/qfxBTAzccouxPYOODidHj9Rx8IDDE4IcNc0AJCVHMTY1gWuuG8uY1AT1BF1q7jZMSfFuy3606+UgNCyYmRkjmJkxgsaGNjZvOsqG0sOs+vtGgoODSB+fxNQrhjF5yhCiojX3kYj4J92Sf44u1S35TqeLkycaOXjA4fk6criOzk4nZrOJESPjGJuawNjUeMamJuhShfilmuomtmw+xtbNx9i/rxqTyYQ1LYEp04YyZdpQLBaNZRORS0O35AeKzk4aK2uoON7KwcP1HNzv4NBBBy0tHQAkp0QxarSFjFkjGDXawrDhsVqjyt90dhoLv4aHQ0sLRHUNOG5sNPbNg7O94hMiuf5GK9ffaKW+vpVtW4wZ0f/7n+08+8w2Ro2OY+KkIUycnMKIkXG6xCsiPqVQ5A82bSIqM5M19/6SunFTGD3Gwo3zxjFqtIWRo+KIiAjxdQ3lbDZtgsxM45b8++6D0lKjPDPT2M/I8G39/EBMTBhzrx7D3KvH0Nzczo5tlWzfVslbb9gpeXk3UdGhTJyYwoRJyYyfkExklC6zicilpVDkR+ZcNYa5X8r2dTVELrqIiBAyZo0kY9ZIOjudHNhXw86umdRL1x/CZIIxY+OZMDGFceOTGDkqDrNZc2eJyMWlUORHNGGiDEZmc5BnDqTbPjoJh6PZs9TMG6+V8/KLuwgLC8aalsC48Umkj09iyNAY3UggIhecQpGI+BWLJYI5V41hzlVj6Ox0cuhgrTED++6TvPDfnXR2OomOCWNceiLj0pNITUsgKUl3W4rI+VMoEhG/ZTYHeeZCmpedTntbJ/sq3MvUnGTThiO4XBAdE0Zqajyp1gTGWhMZPiJWg7ZF5JwpFPmDadNY/tV/Mid1vK9rIgM1bRocP27caXb8OFgsRnnPfTlvIaFm0ickkz7BWEuuubmd/ftqqLBXU1FezfNdPUlhYcGMGesOSQmMGm0hNHRw3gEoIv2nUOQPQkJoirLgCtZdZgErJKR70dfIyO7yAF0INlBERIQwcVIKEyelAMZaf4cO1lJhN5a5cY9JCgoyMWx4LKPHWBg12sLoMfEkJUepN0lEvCgU+YPycu4p+j41E38CjPN1bWQgysvh61+Hr30NfvUr+OUvjfKvf93YT0vzZe0GjZAQM6nWBFKtCYAxCeqxo/Xs31fDwQM1lO+t4t21+wEIDw/2BKTRYyyMGmMhOlqToIoMZgpF/qC2lgl73md9Y4OvayIDVVsL//0vLFxobB9+2CjvuS+XXFCQieEjYhk+IparrhkDGJfcDh1wcGC/gwP7a3jv3f2sKdkDQEJCJCNGxTFiZBwjRsYyYmScgpLIIKJQJCKDSkREiNe4JJfLRU11Mwf213DwgIPDh+t4fc1ez4zycZZwIySN6A5LsXHhuttN5DKkUCQig5rJZCIhMZKExEhmZIwAjMtuNdVNHD5U2/VVxztv76OxsQ2A6OhQRoyMY/iIWIYOi2XosBiSU6I115hIgFMoEhE5RVCQicSkKBKTorhixnDA6FGqdbT0CEq1bCg9jMNRDoDZbCI5JZqhw2IYNtwISsOGxRJnUa+SSKBQKPIHI0bw0rxFxCQN9XVNZKBGjICf/xymTze2I4weB699CWgmkwlLfASW+AimTOv+t9rU1Ebl0XqOHq3n2NF6jh2tY+f2457LbxERIQwZGsOw4TEMHRbLkKHRpAyJ1lglET+kUOQPhgzh3Tl3c1NCkq9rIgM1ZAh84xvG/tSp3eXuMrlsRUaGepYpcXO5XDhqmj1B6eiROuzl1bz/7gGcTpfxvKgQUlJiSBlihKQhXds4S4SmChDxEYUif1BTw+QdbxJ8pea0CVg1NWCzwaxZsH49ZGUZ5TabsR8f79v6ySVlMpmIT4gkPiGSyVOGeMo7OpxUnWzkeGUDxysbqKys59BBBxtKD9He7gSMCSpTUqI9Ycm9n5gUSXCwJqAUuZgUivxBRQUL//0I62+5EsjwdW1kICoqYMECeOopuO8+KC01yhcsMPYVigRj0echQ2MYMjTGq9zpNHqWKivrPYHpeGUDO3ccp7mpHQCTCSzxESQnR5OUHElSUhSJyVEkJ0cTnxCB2axB3iLnS6FIRMTHgoK674CbNLm7Z8nlctHQ0MbxynpOnmg0vk42Ur63mg/eO0hHh9Pr+UlJUSQle39ZdDlOpN8UikRE/JTJZCImJoyYmDDSxnmPOXQ6XdQ6mj1ByR2adu44TvXbTZ6xS8HBQSQkRHpCl/srMdG4vBcerj8DIm761yAiEoCCgrrHLbknonTr7HTiqDEC04kTjVRXNVFV1YS9vIp1Hxykva3Tc2xUVKh3UOraJiRGEhcXrstyMqgoFPmDiAiODEnDGRbu65rIQEVEwMyZEBdnbCMijPKe+yKXiNkc5JlnacIk78fcl+Sqq5o8X1VVTVRXNbK/opra2hZcRicTQUHGNARGr1IE8fGRnmkJ4uMjiLMoNMnlRaHIH0yaROEDv+emMVoMNmBNmgRlZcb+Rz/aXe4uE/ETPS/JjRnb+waAjo5Oamqae4Wmw4fq2Lal0jOrt/FaEBsXjsUSQXxCRPe2KzRZ4nV5TgKLzlYREfEIDjaTnBxNcnJ0n4+3tXXiqGnGUdNMTU0TNZ79ZvZX1FBb2+IZzwTG5JWWeO/QFBcXTpzF2MbGhWmqAfEbCkX+YMMGvrv8djaNLYab03xdGxmIDRtg7lz4y1/gM5+B994zyufONfZnzvRl7UQumNBQs2cOpb44nS7qalu8wpI7QJXvOUmNo5m21k6v50RHhxohyRJuBCZ3aOrxfWiY/lzJxaezzB+4XAR3tuO5kC+Bx+WCtjZwOo2tuy177osMAu5xSJb404+la2lpp9bRYnzVNuNwtFBXa3y/r6KaWkcLTV3zM7lFRIR4QlJsXLixb4nA0vV9TGwYkZGhmn5AzotCkYiIXFLh4SGEDw3pNYllT+1tndTWGqHJCE/dIerokTp27jhOQ32r1/85goJMxMSGERsbTmxsGDGxRliKjQsnNsb4PjY2jOiYMIUn6ZNCkYiI+J2QULNnAsrT6ex0UlfbQl1dK/V1xrauroX6ru2B/TXU1bbQ2NjmFZ5MJoiOMcKTEaLcgckdprof0911g4tCkYiIBCSzOcgzV9OZdHY6aahv7QpPRmDqDk+tHDlcR92O4zTUt3kNEgdj4d7oaKN3KSYmjOjoUGPr/urxfUiIBowHOoUifzBpEr978E9M1y35gWvSJNi6FYYNM7ZWq1Hec19EfMJsDuoauH3mOcOcTheNjW2eXqf6uhbq69toqG+lvr6VhvpWjh2po76hlabG9l7PDw8P7gpKYcTEhHr2jUAV2h2sYsIIDTVjMukSnr9RKPIHERGcSB6ryRsDWUQETJli7CckdJe7y0TE7wUFdc/hNHzEmY/t7HTS2NBGQ4M7MLUZ2wYjPDXUt3HyZA0N9a00NvbugQoJCfKEpKioUKKiw4iKCiEq2uh9MspCPY+FhSlEXQoKRf5g/37mP/9L2qd9G9At+QFp/3740Y/g/vvhiSfgu981yn/0I2N/zBjf1k9ELiizOcgYwB139v/MOp0umpraaKjvGaK6vhraPIv+Nja20djQRktLR6/XCA4OOiUoGdvo6DBjv0d5dFQYEZEhGkw+AAEZiux2O/n5+VitVhwOBxaLhYKCgn4912azUVRUhMViwW63k52dzeLFiy9yjc+iqorMzS9TWvd539ZDBq6qCh5/HG66ydh+4QtGuXtfoUhk0AoKMhmX0aLDgNPfcefW0eGksbGNpkYjMDU2tNLY0GaEpq6yutoWjh6pNx47ZSA5GIPJIyNDPYEpOiqUyK7QFBkZSmRkiPF9VPd+RISCVMCFIrvdTmZmJitXriQnJweA3NxccnNzKSoqOuNzi4uLWbRoETU1NQA4HA4yMzMpLy/vd6gSERG5mIKDgzyTVvaH0+miubnd09PU2GD0QHlCVVd5VZWDpqZ2mpraek2gCUaQCo8IITIypCssdQemyKgQ43t3iIoMJSrKeOxyGh8VcKEoPz+fhIQETyACKCgoIC0tjeLiYq/yvp67YMECz/cWi4W8vDzy8/PJy8vDqgGxIiISYIKCTMals6hQSOnfczo6OmlqNAJSY9e2qamdpq4eKvd+dXUzhw/V0tj1/aljo8C4lNgzPEVF9ghRnnKjJyoyMoSIrlAVEhLkd2EqoEKRw+GguLi41+Uuq9WKxWKhsLDwtKGorKzM08vUU05ODvn5+RQWFqq3SEREBoXgYDOxceZ+jYlyc7lctLZ20NRo9Eo1N7XT2NQVojwByyg/eqTZCFuNfY+RMuoQ5BWUIrpC1B0fm+qzhYQDKhTZbDYA0tJ6D0a2Wq2ex/tit9v7LE/oulOozJermQ8ZwptzFxASn+S7Osj5GTIEHnoIJkwwtkOGGOU990VEApjJZDJmIw8PISHxzHND9dTZ6aS5uZ3mpnaamtppbjYCVHNjO01d5c1dl/VOnmjEbPZd71FAhSJ3sLFYLL0ec4cb98Dr0ykvL/f63n3s6ULTqZxOJ3V1dV5lYWFhhIWF9ev5fRoxAtuN9zMveejAX0N8a8QIWLbM2J81q7vcXSYiMkiZzUE9Bpr7t4Ccvzyh5zwwXc4WbjIyMgBO25t0piDV0549e4iLi/P6Wna+f/jq6xm7fzPmpobzex3xnfp6eP11OHLE2NbXd5fV1/u4ciIi0h8BFYrcA6Grq6t7PeZwOIDThxur1crixYspKytj+fLlnufk5+cDfQetvqSnp1NbW+v1tXTp0nP8SU6xZw8P/COfiEP7zu91xHf27DFux3/tNWO7Z0932Z49vq6diIj0Q0BdPnOHIncA6skdlM50B1lhYSFpaWmsWrWKwsJCsrKyPAOvc3Nz+1WHoKAgYmNjz7HmIiIi4u8CMhSdOi4IjMtm7ktkZ7JkyRKWLFni+T47OxvA61Z9ERERGXwC6vKZxWIhJyen17ggu92Ow+HodRnrbIOnbTYbNpuNwsLCfo8pEhERkctTQIUiMCZqrK6u9gpG+fn5ZGVlec1RtHz5cs+Ejn2x2+3k5uayePFi3y/zERJCbUwiruAQ39ZDBi4kxLgDLTzc2IaEdJeFqF1FRAJBQF0+A+MSWmlpKfn5+ZSUlGC327Farb0mXnRP6NjXGCP3ch9Lly71upTmM9Om8fMvPsU8a7qvayIDNW0aHDpk7N99d3e5u0xERPxewIUiMALP2dY5y8nJ6TW79YoVKygtLSUtLY2Kigq/u2Tmb9Odi4iIDCYBGYoGyueXyU5nyxa++fv72JP+JNyo9dcC0pYtcNtt8Otfw1e/Ci++aJTfdpuxP22ab+snIiJnNahCkd9qbyeuvoqgjnZf10QGqr0dDh+GlhZj297Vlj33RUTErwXcQGsRERGRi0GhSERERASFIhERERFAY4r8Q3o6j99TwMSRY31dExmo9HRj3bPx441tetf0Cj33RUTErykU+YOYGPaNvoL0yGhf10QGKiYGbrzR2B8+vLvcXSYiIn5Pl8/8weHDZL/xBKEnjvm6JjJQhw/D0qWwfr2xPXy4u+zwYV/XTkRE+kGhyB9UVnL9e6sJrTnp65rIQFVWwqOPwq5dxraysrusstLXtRMRkX5QKBIRERFBoUhEREQEUCgSERERAXT3mX9ITGT9FbfSGRvv65rIQCUmwmc/C2PHGtvERKO8576IiPg1hSJ/MGYM/7ntq9wydISvayIDNWYMPPaYsX/NNd3l7jIREfF7unzmD5qbSTmxn6DWFl/XRAaquRm2bYPqamPb3Nxd1tzs69qJiEg/KBT5gx07+MqfP0/k/r2+rokM1I4dMHUqvPiisd2xo7tsxw5f105ERPpBoUhEREQEhSIRERERQKFIREREBNDdZ/7BZKLDHIzLZPJ1TWSgTCYIDYWgIGPrbsue+yIi4tcUivzBzJk8/H/Pcev4Cb6uiQzUzJnQ2mrs33NPd7m7TERE/J4un4mIiIigUOQfduzgC3/5MhG6JT9w7dgBGRnwv/8ZW/ct+e59ERHxewpF/qC5meGV5Zq8MZA1N8OGDVBba2zdkze690VExO8pFImIiIigUCQiIiICKBSJiIiIALol3z+kpvKP+UsZO2yUr2siA5WaCqtXw6xZxjY11SjvuS8iIn5NocgfxMezbeJ1jIyJ83VNZKDi4yE319jvGYLcZSIi4vd0+cwfVFZy9bpnCKk+6euayEBVVsIvfgFbtxrbysrusspKX9dORET6QaHIHxw+zIdfe4ywk8d8XRMZqMOH4ZvfhE2bjO3hw91lhw/7unYiItIPCkUiIiIiKBSJiIiIAApFIiIiIoDuPvMPcXHsSJtDZ1SMr2siAxUXB7ffDsOGGdu4rjsJe+6LiIhfUyjyB2lpPHX397ltxBhf10QGKi0NnnvO2L/55u5yd5mIiPg9XT7zB+3tRDbVYupo93VNZKDa2+HECWhqMrbt7d1l7WpXEZFAoFDkD7Zs4Tu/u4dI+y5f10QGassWSEmBf//b2G7Z0l22ZYuvayciIv2gUCQiIiKCQpGIiIgIoFDkX0wmX9dARERk0FIoEhEREUG35PuH6dP54VeKuDltkq9rIgM1fTrU1kJ4uDE3UVSUUV5b270vIiJ+TaHIH5jNtIZFgtns65rIQJnNEBtr7IeGdpe7y0RExO/p8pk/2LOHzxR9l/BDFb6uiQzUnj1w661QUmJs9+zpLtuzx9e1ExGRflAo8gf19aTvK8Pc1OjrmshA1dfDK6/A8ePGtr6+u6y+3te1ExGRflAoEhEREUGhSERERARQKBIREREBdPeZfxg1iufmfZ6UlGG+rokM1KhR8LvfQUaGsR01yijvuS8iIn5NocgfJCfzfsZHuc2S6OuayEAlJ8MXv2jsT+ox35S7TERE/J4un/mD6mqmb3uVkDqHr2siA1VdDU89BeXlxra6urusutrXtRMRkX5QKPIH+/ax4IWfE3bskK9rIgO1bx988pPw3nvGdt++7rJ9+3xbNxER6ReFIhEREREUikREREQAhSIRERERQHef+YeoKA4Mm4gzPNLXNZGBioqCuXMhIcHYRkUZ5T33RUTErykU+YMJE/jTfT/nI6Otvq6JDNSECfDuu8b+bbd1l7vLRETE7+nymYiIiAgKRf6hrIxHfvoRonZv9XVNZKDKysBkgqefNrZlZd1lZWW+rp2IiPSDQpGIiIgICkUiIiIigEKRiIiICKBQJCIiIgLolnz/MHkyP/vsSq4em+7rmshATZ4Me/ZAUpKxHTnSKO+5LyIifi0gQ5Hdbic/Px+r1YrD4cBisVBQUNCv59psNgoLCz3Pra6uZuXKlVgslotb6TMJD6c6fjiusDDf1UHOT3g4jBtn7Pc8l9xlIiLi9wLu8pndbiczM5OFCxdSUFBAYWEhdrud3Nzcsz7XZrORm5tLUVGR57mzZ88mMzPzEtT89Fp37iTnf8sJPljh03oEitbWVh5++GFaW1t9XZVuFRVw333w5pvGtqKiu6zCd+3ql5+Vn9Jn1X/6rPpPn1X/+cVn5QowOTk5LqvV6lVWXl7uAlxFRUVnfG5WVpYrJyenz+eWl5f36/1vv/32c6twP9S/8YbLBa73f198wV/7clRbW+sCXLW1tb6uSrfSUpcLXK6nnjK2paXdZaWlPquWX35WfkqfVf/ps+o/fVb9d7E/q/78/Q6oniKHw0FxcTFZWVle5VarFYvFQmFh4Vlfw2azeX1f1jWxXkJCwoWrqIiIiAScgApF7kCTlpbW6zGr1dor8JwqLy8Ph8NBZmYmDocDh8PBsmXLKCoq8u2YIhEREfG5gBpobbfbAfoMMO6eHvfA677k5ORQWFhIXl4eqampzJo1i6KiIqzW/i/E6nQ6qaur8yoLCwsjTIOkRUREAprJ5XK5fF2J/lq+fDn5+fkUFRWRk5Pj9Vhubi7FxcWUlpaSkZFxxtdxHwtQUFDAkiVL+l2H5ORkTp486VU2fvx4JkyY0O/XOJXT6WTPnj2kp6cTFBRQnXc+oc+r//RZ9Z8+q/7TZ9V/+qz672J/VuXl5Wzbtu2MxwRUT5G7R6e6urrXYw6HA+i7F6kn96385eXlZGdnk5+fT3l5eb/GIwEcOnSo18h49RSJiIgEvoAMRe4A1JM7KJ3pUlh+fj7FxcWUl5cDUFpayrx581ixYgW5ubm9BnD3RQFIRETk8hRQfXnuwOMONT3Z7fazXjZbsWKF1zEWi4U1a9YAUFJScgFrKiIiIoEmoEKRxWIhJyen111mdrsdh8PB0qVLe5X35J7F+tTXtFgsJCYmXpQ6i4iISGAIqFAExsDo6upqr2CUn59PVlaW1+Dr5cuXk5aW5hlQDbB06VJsNptXWHLPU7R48eJLUHsRERHxVwE1pgiM3p7S0lLy8/MpKSnBbrdjtVp7rX3mntCx5xijnJwcioqKyMvL8zxut9spLS3VPEUiIiKDXEDdki8iIiJysQTc5TORU506dkzkQtM5JheTzi//oVDkY3a7ndzcXPLz88nLyyM/P9/XVfJraWlpmEwmr6/BzOFwkJ+fT3Z29mmP0Tlm6M9nBTrH3Gw2G5mZmZhMJtLS0li+fHmfx+n86v9nBTq//P28CrgxRZcTu91OZmYmK1eu9AwSz83NJTc3l6KiIh/Xzv8UFxdjtVq9BtSnpaWd0zItlxObzUZJSQnLly8/7XQUOscM/fmsQOeYW3FxMcuWLWPhwoXMmjWL1atXk5+fT1VVldf4TZ1f/f+s3McO5vPLZrNRWFjo+VwKCgrIz8/HYrF43ezk0/PKJT6Tk5PjslqtXmXl5eUuwFVUVOSjWvmvrKwsV01Nja+r4XcsFosrIyOjz8d0jnk702flcukcczv1M6qpqXFZLBbXqX8ydH71/7NyuXR+LVmypFcZ4MrJyfEq8+V5pVDkIzU1NS7AtXjx4l6PWSwWV1ZWlg9q5b9KS0s9n1dRUdGg/sVyKqvV2ucfep1jvZ3us3K5dI65lZaWukpKSnqVL1682AW4ysvLXS6Xzi+Xq/+flftYnV+9WSwWV2Fhoed7X59XGlPkI+55ltLS0no9ZrVae01QOditWrUKwLMkS3x8PCtWrPBxrfybzrFzo3PMkJGR0eeSR+5pS9yXenR+9f+zAp1ffbHZbGRlZXldOvP1eaVQ5CPuuw36mh8pISEB6HuNt8GqoKAAl8tFeXk5S5YsASAvL2/Q/1I5E51j50bn2JmVlZV5jYXR+XV6p35WoPPrVCtWrCA7O5vZs2d7lfv6vFIo8jF3I/fkPhl0m2Zv7ok63RNu5uXl+bpKfk/n2LnROdab3W7HbrezcuXKXo/p/PJ2ps8KdH6BseJESUkJFouF/Px8cnNzex3jq/NKochH3N2q1dXVvR5zp2DNsn16GRkZnjsYBuMv3v7QOXZ+dI51y8vLo6ioyOt80fnVt74+q74M5vNryZIlFBUVUVNTQ0ZGBsXFxZ7PwNfnlUKRj7gbvq9uQPfJMFhu0xwo97X8wdpFfzY6x86fzjE8/5M/dSoDnV+9ne6zOh2dX3h61E4NRb46rxSKfMTdqOXl5b0es9vt/f5HJeizOg2dYxfOYP2siouLSUtL63PBbJ1f3s70WZ3NYPusenKfR6dufXVeKRT5iMViIScnp9dIervdjsPhYOnSpT6qWeAoLi7uNZhRuukcO3+D+Ryz2WxUV1f3+iPvPp90fnU722d1OoP5/HJbv349VqvVE4Z8fl5d1Bv+5YzKy8tdFovFa56LnJycQTG/x7koLS11Wa1WV0FBgWduj9LS0kE/EZrbmSYk1Dnm7XSflc4xbyUlJZ7Po+fX4sWLvSba0/nVv89K55chKyur15xEWVlZXvM5uVy+Pa9MLpfLdXFjl5yJ3W4nPz8fq9WK3W733Jkg3RwOB7m5uaxfvx6AxYsXD7ib+nJis9koKiry3NJbUFBAVlZWr+5lnWNn/6x0jnUrKysjMzPztI8XFRX1ujV/sJ5f/f2sdH4ZcnNzsdlsJCQkkJOTQ2Jiomd6glP56rxSKBIRERFBY4pEREREAIUiEREREUChSERERARQKBIREREBFIpEREREAIUiEREREUChSERERARQKBIREREBFIpEREREAIUiEREREUChSERERARQKBIREREBFIpERLzk5eX5ugoi4iPBvq6AiIi/yM3NxW63+7oaIuIj6ikSkUHPbreTmZlJcXExWVlZvq6OiPiIQpGIDGo2m43c3FzKysoAKC4uJjMz0/O9iAweJpfL5fJ1JUREfC0+Ph6Hw4F+JYoMXuopEpFBz26343A4dOlMZJBTKBKRQc9mswGQnZ3t45qIiC8pFInIoFdSUgJATk6Oj2siIr6kMUUiMujFx8cDUFNT4+OaiIgvqadIRAa1srIyjScSEUChSEQGOY0nEhE3hSIRGdTc44mysrKw2Wyan0hkEFMoEpFBzd1TZLPZcDgcZGRk+LhGIuIrCkUiMqhlZWVhsVhwOBy6+0xkkNPdZyIiIiKop0hEREQEUCgSERERARSKRERERACFIhERERFAoUhEREQEUCgSERERARSKRERERACFIhERERFAoUhEREQEUCgSERERARSKRERERACFIhERERFAoUhEREQEgP8P/LIt45ad8ZkAAAAASUVORK5CYII=",
"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",
"plt.plot(t_approx, Ft_approx, label='approx')\n",
"images_rel_td = lens.img_td - min(lens.img_td)\n",
"for imtd in images_rel_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.ylim(bottom=0.8)\n",
"# plt.xlim(6, 8)\n",
"# plt.xscale('log')\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9a4922b1-e8b5-435b-9dad-ca3ec80b1692",
"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
}
notebooks/Fw_computation.ipynb 0 → 100644
This source diff could not be displayed because it is too large. You can view the blob instead.
utils/generate_ampfac_charged_lens.py 0 → 100644
#!/home/uddeepta.deka/.conda/envs/igwn-py38-lw/bin/python3
__author__ = "U.Deka"
"""
This code is used to generate amplification factor in the time
and frequency domain for the charged lens model using
- contour integration method
- area method using histogram
"""
import sys
sys.path.append("../utils/")
import numpy as np
import matplotlib.pyplot as plt
import warnings
from lenses import ChargedLens
import time
import os
warnings.filterwarnings("ignore")
#===========================================================================
#============================ Histogram method =============================
#===========================================================================
def Ftilde_hist(lens, t_max, dt, dx_fac=1, verbose=False, dx1=None, dx2=None):
"""
Computes F(t) with the histogram method.
Params:
---------------
lens: object of a charged lens class
t_max: maximum time delay to be resolved (relative to global minima)
dt: bin width of histogram
"""
# We assume that atleast one image exists
assert len(lens.img)!=0
minima_td = min(lens.img_td)
t_max += minima_td
x_max = 1.2 * np.sqrt(2 * t_max)
dtdx1, dtdx2 = lens.gradTimeDelay(x_max, x_max)
if not dx1:
dx1 = dt / dtdx1 * dx_fac
if not dx2:
dx2 = dt / dtdx2 * dx_fac
x1_ = np.arange(-x_max + lens.y1, x_max + lens.y1, dx1)
x2_ = np.arange(-x_max + lens.y2, x_max + lens.y2, dx2)
X1, X2 = np.meshgrid(x1_, x2_)
td_samples = lens.timeDelay(X1, X2).ravel()
td_samples = td_samples[td_samples >= minima_td]
td_samples -= np.min(td_samples)
bins = int(max(td_samples) / dt)
td_hist, bin_edges = np.histogram(td_samples, bins=bins)
t_mid = np.zeros_like(bin_edges)
t_mid[1:] = (bin_edges[1:] + bin_edges[:-1]) / 2.0
Ft = np.zeros_like(bin_edges)
Ft[1:] = dx1 * dx2 * td_hist / dt / 2.0 / np.pi
Ft[0] = np.sqrt(lens.img_mag[np.argmin(lens.img_td)])
if verbose:
print(f"Tmax = {t_max} | dt = {dt} | xmax = {x_max} | dx1 = {dx1} | dx2 = {dx2}")
print(f"Number of td samples = {len(td_samples)}")
cut_idxs = t_mid<=t_max
return t_mid[cut_idxs], Ft[cut_idxs]
#===========================================================================
#============================== Contour method =============================
#===========================================================================
def Ftilde_contour(lens, t_max, dt, dx_fac=1, verbose=False, dx1=None, dx2=None, t_min=0):
"""
Computing the raw F(t)
for charged lens using contour summation method
Params:
---------------
lens: object of a charged lens class
t_max: maximum time delay to be resolved (relative to global minima)
dt: bin width of histogram
"""
# We assume that atleast one image exists
assert len(lens.img)!=0
minima_td = min(lens.img_td)
t_max += minima_td
x_max = 1.2 * np.sqrt(2 * t_max)
dtdx1, dtdx2 = lens.gradTimeDelay(x_max, x_max)
if not dx1:
dx1 = dt / dtdx1 * dx_fac
if not dx2:
dx2 = dt / dtdx2 * dx_fac
x1_ = np.arange(-x_max + lens.y1, x_max + lens.y1, dx1)
x2_ = np.arange(-x_max + lens.y2, x_max + lens.y2, dx2)
X1, X2 = np.meshgrid(x1_, x2_)
n = int((t_max-minima_td) / dt)
Tarray = np.logspace(0, np.log10(t_max+1), n)
Tarray -=1
Tarray = Tarray + t_min + minima_td
plt.figure()
contour_levels = plt.contour(X1, X2, lens.timeDelay(X1, X2), Tarray)
level_coords = contour_levels.allsegs
Ft = []
for level in level_coords:
res_sum = 0
for i in range(len(level)):
if len(level[i])!=0:
x_coord, y_coord = np.vstack(level[i]).T
ds = np.sqrt(np.diff(x_coord)**2 + np.diff(y_coord)**2)
gtd = lens.gradTimeDelay(x_coord, y_coord)
gradnorm = np.linalg.norm([gtd[0],gtd[1]], axis = 0)[:-1]
res = np.divide(ds, gradnorm)
res_sum += np.sum(res) / 2. / np.pi
Ft.append(res_sum)
plt.close()
Tarray-=minima_td
Ft[0] = np.sqrt(lens.img_mag[np.argmin(lens.img_td)])
if verbose:
print(f"x_max = {x_max} | t_max = {t_max} | dx1 = {dx1}, dx2 = {dx2}")
print(f"Number of td samples = {len(Tarray)}")
return Tarray, np.array(Ft)
#===========================================================================
#=========================== Peak reconstruction ===========================
#===========================================================================
class ReconstructPeak:
"""
Class for reconstructing peak of F(t).
"""
def __init__(self, lens, t, Ft, peak_window_frac=0.05):
self.lens = lens
self.t = t
self.Ft = Ft
self.peak_window_frac = peak_window_frac
def get_approx_Ft(self):
assert len(self.lens.img)!=0
image_timedelays = self.lens.img_td - min(self.lens.img_td)
Ft = np.zeros_like(self.t)
for ii in range(len(self.lens.img_td)):
img_morse_idx = self.lens.img_midx[ii]
if img_morse_idx == 0.5:
Ft -= (np.sqrt(np.abs(self.lens.img_mag[ii])) *
np.log(np.abs(image_timedelays[ii]-self.t)) / np.pi)
elif (img_morse_idx==0) or (img_morse_idx==1):
sel_idxs = self.t > image_timedelays[ii]
Ft[sel_idxs] += (-1)**(img_morse_idx) * np.sqrt(self.lens.img_mag[ii])
return np.asarray(Ft)
def Tj(self):
peak_idx = np.argmax(-1 * np.asarray(self.lens.img_mag))
return self.lens.img_td[peak_idx] - np.min(self.lens.img_td)
def find_peak_window(self):
t_peak = self.Tj()
t_start = (1- self.peak_window_frac) * t_peak
t_end = (1 + self.peak_window_frac) * t_peak
start_of_window = np.where(self.t >= t_start)[0][0]
end_of_window = np.where(self.t >= t_end)[0][0]
real_peak_loc = np.where(self.t >= t_peak)[0][0]
return [start_of_window, end_of_window, real_peak_loc]
def get_peak_reconstructed_Ft(self):
left_idx, right_idx, estimated_real_peak = self.find_peak_window()
# get approx peak correction from analytical expression
Ft_approx = self.get_approx_Ft()
# Get the difference between the Ft and Ft_approx on the left and right
diff_left = self.Ft[left_idx] - Ft_approx[left_idx]
diff_right = self.Ft[right_idx] - Ft_approx[right_idx]
# peak of the approximated Ft
peak_Ft_approx = np.argmax(Ft_approx)
# Rescale the left and right part
Ft_approx_rescaled_left = Ft_approx[:peak_Ft_approx] + diff_left
Ft_approx_rescaled_right = Ft_approx[peak_Ft_approx:] + diff_right
# combine to get peak corrected Ft
Ft_corrected = np.copy(self.Ft)
Ft_corrected[left_idx : peak_Ft_approx] = Ft_approx_rescaled_left[left_idx : peak_Ft_approx]
Ft_corrected[peak_Ft_approx : right_idx+1] = Ft_approx_rescaled_right[0 : right_idx + 1 - peak_Ft_approx]
Ft_corrected[peak_Ft_approx] = 0.5 * (Ft_corrected[peak_Ft_approx - 1] + Ft_corrected[peak_Ft_approx + 1])
return Ft_corrected
#===========================================================================
#================== Frequency domain amplification factor ==================
#===========================================================================
def ampFac(t, Ft):
"""
Returns the angular frequency (dimensionless)
and amplification factor
Input params:
---------------
t: time array
Ft: time series data
"""
fft_ = np.fft.rfft(Ft) * (t[1] - t[0])
freq_ = np.fft.rfftfreq(len(t), t[1] - t[0]) * 2.0 * np.pi
ampfac = fft_ * freq_ * 1j
ampfac = np.conj(ampfac) * np.exp(1j * freq_ * t[0])
ampfac += Ft[-1]
return freq_, ampfac
def Fw_geom(lens, w):
"""
Returns the amplification factor in
geometric optics limit.
Input params:
---------------
lens: object of class ChargedLens
w : list of angular frequency
"""
assert len(lens.img)!=0
delta_img_td = lens.img_td - min(lens.img_td)
F_GO_list = []
for ii in range(len(delta_img_td)):
F_GO_list.append(np.sqrt(abs(lens.img_mag[ii]))*np.exp(1j*w*delta_img_td[ii] -
1j*np.pi*lens.img_midx[ii]))
return np.sum(F_GO_list, axis=0)
#===========================================================================
#=============================== Driver code ===============================
#===========================================================================
def usage():
print(
"Computes the lensing amplification factor in time and frequency domain due to charged lens."
)
print("")
print("The following options can be used to overwrite default parameters")
print("\t-y: Impact parameter, assuming y2=0 [Default: 0.1]")
print("\t-q: Effective charge [Default: 0]")
print("\t-t_max: max absolute time to integrate [Default: 30]")
print("\t-verbose: flag to print output")
print("\t-dt_fac: Factor to downscale time spacing [Default: 7.5]")
print("\t-w_min: Minimum frequency to be resolved in freq. domain [default: 0]")
print("\t-w_max: Maximum frequency to be resolved in freq. domain [default: 200]")
print(
"\t-folder_path: Relative path to the folder to store results [default: ampfac_data]"
)
print("\t-method: Method to compute the amplification factor ('contour', 'hist', 'geom') [default:'contour']")
def main():
print(" --------------------------------------------------------")
print(" : generate_ampfac_charged_lens.py :: use flag -h for list of options:")
print(" --------------------------------------------------------")
# set default values for parameters
# impact parameter
y = 0.1
# effective charge
q = 0.
# integrate to time T_MAX
T_MAX = 30
# print output
verbose = False
# dt factor. Larger dt_fac => smaller dt
dt_fac = 7.5
# Minimum dimensionless angular frequency
W_MIN = 0
# Maximum dimensionless angular frquency
W_MAX = 200
# folder path to store data
folder_path = os.getcwd() + "/ampfac_data/"
# computation method to use
method = 'contour'
## now look for flags to overwrite default values
for i in range(len(sys.argv)):
if sys.argv[i] == "-h":
usage()
return
if sys.argv[i] == "-y":
y = float(sys.argv[i+1])
if sys.argv[i] == "-q":
q = float(sys.argv[i+1])
if sys.argv[i] == "-t_max":
T_MAX = float(sys.argv[i+1])
if sys.argv[i] == "-verbose":
verbose = True
if sys.argv[i] == "-dt_fac":
dt_fac = float(sys.argv[i+1])
if sys.argv[i] == "-w_min":
W_MIN = float(sys.argv[i+1])
if sys.argv[i] == "-w_max":
W_MAX = float(sys.argv[i+1])
if sys.argv[i] == "-folder_path":
folder_path = os.getcwd() + "/" + str(sys.argv[i+1]) + "/"
if sys.argv[i] == "-method":
method = str(sys.argv[i+1])
isExist = os.path.exists(folder_path)
if not isExist:
os.makedirs(folder_path)
print("Created new folder to store data.")
## Final params
DT_FAC = 50
DT = np.pi / W_MAX / DT_FAC
## current params
w_max = 20
dt = np.pi / w_max / dt_fac
## determining peak correction window
if y < 0.015:
win_frac = 0.5
elif y < 0.025:
win_frac = 0.4
elif y < 0.1:
win_frac = 0.2
elif y < 0.5:
win_frac = 0.05
elif y < 0.8:
win_frac = 0.02
elif y >=0.8:
win_frac = 0.01
tStart = time.perf_counter()
print(f"Generating data for y = {y}, q={q} using {method} method ... ")
# Creating the lens object
lens = ChargedLens(y1=y, qeff=q)
if method == "contour":
# you may try changing dx1 and dx2 values. But these seem to work well
t_vals, Ft_vals = Ftilde_contour(lens, T_MAX, dt, dx1=0.025, dx2=0.025, verbose=verbose)
# interpolating
num = int(2**np.ceil(np.log2((max(t_vals) - min(t_vals))/DT)))
t_dense = np.linspace(min(t_vals), max(t_vals), num)
Ft_interped = np.interp(t_dense, t_vals, Ft_vals)
# applying peak correction
rp = ReconstructPeak(lens, t_dense, Ft_interped, peak_window_frac=win_frac)
Ft_recons = rp.get_peak_reconstructed_Ft()
# doing the FFT
w_vals, Fw_vals = ampFac(t_dense, Ft_interped)
idxs = (w_vals > W_MIN) & (w_vals <= W_MAX)
w_vals = w_vals[idxs]
Fw_vals = Fw_vals[idxs]
_, Fw_recons = ampFac(t_dense, Ft_recons)
Fw_recons = Fw_recons[idxs]
# saving results to file
Ft_datafile = folder_path + f'FtCont_y{y}_q{q}_wmax{W_MAX}_dtfac{DT_FAC}_tmax{T_MAX}.txt'
Fw_datafile = folder_path + f'FwCont_y{y}_q{q}_wmax{W_MAX}_dtfac{DT_FAC}_tmax{T_MAX}.npz'
header = f"Created using generate_ampfac_charged_lens.py with contour method.\nSettings: t_max={T_MAX}, dt={t_dense[1]-t_dense[0]}\nUsage: \nimport numpy as np\ndata = np.loadtxt('{Ft_datafile}')\nt = data[:,0]\nFt = data[:,1]\nFt_recons = data[:,2]\n"
with open(Ft_datafile, 'w') as f:
np.savetxt(f, np.array([t_dense, Ft_interped, Ft_recons]).T, header=header)
np.savez(Fw_datafile, w=w_vals, Fw=Fw_vals, Fw_reconstructed=Fw_recons)
elif method == "hist":
# you may try changing the dx_factor
t_vals, Ft_vals = Ftilde_hist(lens, T_MAX, dt, dx_fac=1, verbose=verbose)
# interpolating
num = int(2**np.ceil(np.log2((max(t_vals) - min(t_vals))/DT)))
t_dense = np.linspace(min(t_vals), max(t_vals), num)
Ft_interped = np.interp(t_dense, t_vals, Ft_vals)
# applying peak correction
rp = ReconstructPeak(lens, t_dense, Ft_interped, peak_window_frac=win_frac)
Ft_recons = rp.get_peak_reconstructed_Ft()
# doing the FFT
w_vals, Fw_vals = ampFac(t_dense, Ft_interped)
idxs = (w_vals > W_MIN) & (w_vals <= W_MAX)
w_vals = w_vals[idxs]
Fw_vals = Fw_vals[idxs]
_, Fw_recons = ampFac(t_dense, Ft_recons)
Fw_recons = Fw_recons[idxs]
# saving results to file
Ft_datafile = folder_path + f'FtHist_y{y}_q{q}_wmax{W_MAX}_dtfac{DT_FAC}_tmax{T_MAX}.txt'
Fw_datafile = folder_path + f'FwHist_y{y}_q{q}_wmax{W_MAX}_dtfac{DT_FAC}_tmax{T_MAX}.npz'
header = f"Created using generate_ampfac_charged_lens.py with histogram method.\nSettings: t_max={T_MAX}, dt={t_dense[1]-t_dense[0]}\nUsage: \nimport numpy as np\ndata = np.loadtxt('{Ft_datafile}')\nt = data[:,0]\nFt = data[:,1]\nFt_recons=data[:,2]\n"
with open(Ft_datafile, 'w') as f:
np.savetxt(f, np.array([t_dense, Ft_interped, Ft_recons]).T, header=header)
np.savez(Fw_datafile, w=w_vals, Fw=Fw_vals, Fw_reconstructed=Fw_recons)
elif method == "geom":
num_w = 2048
w_arr = np.logspace(np.log10(1e-1), np.log10(W_MAX), num_w)
Fw_g = Fw_geom(lens, w_arr)
Fw_datafile = folder_path + f'FwGeom_y{y}_q{q}_wmax{W_MAX}_num{num_w}.npz'
np.savez(Fw_datafile, w=w_arr, Fw=Fw_g)
else:
print("Unknown method")
return
tEnd = time.perf_counter()
print(f"Data saved in {folder_path}\nTime taken to generate ampfac:{tEnd - tStart}s.")
if __name__ == "__main__":
main()
\ No newline at end of file
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or sign in to comment