Notebooks/Simulations/simulating_different_telescopes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulating different telescopes\n",
    "This notebooks provides examples in how to use the lenstronomy.SimulationAPI modules in simulating (realistic) mock lenses taylored to a specific observation and instrument and makes a montage of different telescope settings currently available.\n",
    "\n",
    "The module enables to use the astronomical magnitude conventions and can translate those into the lenstronomy core module configurations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import copy\n",
    "import os\n",
    "import numpy as np\n",
    "import scipy\n",
    "from PIL import Image\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "# make sure lenstronomy is installed, otherwise install the latest pip version\n",
    "try:\n",
    "    import lenstronomy\n",
    "except:\n",
    "    !pip install lenstronomy\n",
    "\n",
    "# lenstronomy module import\n",
    "from lenstronomy.Util import image_util, data_util, util\n",
    "import lenstronomy.Plots.plot_util as plot_util\n",
    "from lenstronomy.SimulationAPI.sim_api import SimAPI\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define camera and observations\n",
    "As an example, we define the camera and observational settings of a LSST-like observation. We define one camera setting and three different observations corresponding to g,r,i imaging.\n",
    "\n",
    "For the complete list of possible settings, we refer to the SimulationAPI.observation_api classes. There are pre-configured settings which approximately mimic observations from current and future instruments. Be careful using those and check whether they are sufficiently accurate for your specific science case!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Instrument setting from pre-defined configurations\n",
    "\n",
    "from lenstronomy.SimulationAPI.ObservationConfig.DES import DES\n",
    "from lenstronomy.SimulationAPI.ObservationConfig.LSST import LSST\n",
    "from lenstronomy.SimulationAPI.ObservationConfig.Euclid import Euclid\n",
    "from lenstronomy.SimulationAPI.ObservationConfig.Roman import Roman\n",
    "\n",
    "DES_g = DES(band='g', psf_type='GAUSSIAN', coadd_years=3)\n",
    "DES_r = DES(band='r', psf_type='GAUSSIAN', coadd_years=3)\n",
    "DES_i = DES(band='i', psf_type='GAUSSIAN', coadd_years=3)\n",
    "des = [DES_g, DES_r, DES_i]\n",
    "\n",
    "LSST_g = LSST(band='g', psf_type='GAUSSIAN', coadd_years=10)\n",
    "LSST_r = LSST(band='r', psf_type='GAUSSIAN', coadd_years=10)\n",
    "LSST_i = LSST(band='i', psf_type='GAUSSIAN', coadd_years=10)\n",
    "lsst = [LSST_g, LSST_r, LSST_i]\n",
    "\n",
    "Roman_g = Roman(band='F062', psf_type='PIXEL', survey_mode='wide_area')\n",
    "Roman_r = Roman(band='F106', psf_type='PIXEL', survey_mode='wide_area')\n",
    "Roman_i = Roman(band='F184', psf_type='PIXEL', survey_mode='wide_area')\n",
    "roman = [Roman_g, Roman_r, Roman_i]\n",
    "\n",
    "# lenstronomy provides these setting to be imported with the SimulationAPI.ObservationConfig routines.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import real galaxy to be the source galaxy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1788321.7625399919\n"
     ]
    }
   ],
   "source": [
    "# import NGC1300 jpg image and decompose it\n",
    "import imageio.v2 as imageio\n",
    "from scipy.ndimage import gaussian_filter\n",
    "\n",
    "# find path to data\n",
    "path = os.getcwd()\n",
    "dirpath, _ = os.path.split(path)\n",
    "module_path, _ = os.path.split(dirpath)\n",
    "ngc_filename = os.path.join(module_path, 'Data/Galaxies/ngc1300.jpg')\n",
    "# read data, this only works if you execute this notebook within the environment of the github repo!\n",
    "ngc_data = imageio.imread(ngc_filename, mode='F', pilmode=None)\n",
    "\n",
    "# subtract the median of an edge of the image\n",
    "median = np.median(ngc_data[:200, :200])\n",
    "ngc_data -= median\n",
    "\n",
    "# resize the image to square size (add zeros at the edges of the non-square bits of the image)\n",
    "nx, ny = np.shape(ngc_data)\n",
    "n_min = min(nx, ny)\n",
    "n_max = max(nx, ny)\n",
    "ngc_square = np.zeros((n_max, n_max))\n",
    "x_start = int((n_max - nx)/2.)\n",
    "y_start = int((n_max - ny)/2.)\n",
    "ngc_square[x_start:x_start+nx, y_start:y_start+ny] = ngc_data\n",
    "\n",
    "# we slightly convolve the image with a Gaussian convolution kernel of a few pixels (optional)\n",
    "sigma = 5\n",
    "ngc_conv = gaussian_filter(ngc_square, sigma, mode='nearest', truncate=6)\n",
    "\n",
    "# we now degrate the pixel resoluton by a factor.\n",
    "# This reduces the data volume and increases the spead of the Shapelet decomposition\n",
    "factor = 25  # lower resolution of image with a given factor\n",
    "numPix_large = int(len(ngc_conv)/factor)\n",
    "n_new = int((numPix_large-1)*factor)\n",
    "ngc_cut = ngc_conv[0:n_new,0:n_new]\n",
    "x, y = util.make_grid(numPix=numPix_large-1, deltapix=1)  # make a coordinate grid\n",
    "ngc_data_resized = image_util.re_size(ngc_cut, factor)  # re-size image to lower resolution\n",
    "\n",
    "print(np.sum(ngc_data_resized))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define model settings\n",
    "\n",
    "The model settings are handled by the SimulationAPI.model_api ModelAPI class. \n",
    "The role is to return instances of the lenstronomy LightModel, LensModel, PointSource modules according to the options chosen by the user. Currently, all other model choices are equivalent to the ones provided by LightModel, LensModel, PointSource.\n",
    "The current options of the class instance only describe a subset of possibilities and we refer to the specific class instances for details about all the possibilities.\n",
    "\n",
    "For this example, we chose a single lens plane and a single source plane, elliptical Sersic profile for the deflector, the interpolated Galaxy as the source and an additional lensed point source."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "kwargs_model = {'lens_model_list': ['SIE', 'SHEAR'],  # list of lens models to be used\n",
    "                'lens_light_model_list': ['SERSIC_ELLIPSE'],  # list of unlensed light models to be used\n",
    "                'source_light_model_list': ['INTERPOL'],  # list of extended source models to be used, here we used the interpolated real galaxy\n",
    "                'point_source_model_list': ['SOURCE_POSITION']  # list of point source models to be used\n",
    "    }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Brightness definitions in magnitude space\n",
    "One core feature is the support of light profile amplitudes in astronomical magnitude space (at least for few selected well defined brightness profiles).\n",
    "\n",
    "We first define all parameters in magnitude space and then use the SimAPI routine to translate the arguments into lenstronomy conventions used by the ImSim module. The second model of each light component we defined as 'INTERPOL', which sets an interpolation grid given an image. This can be used to past real galaxies as lenses or sources into lenstronomy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# g-band\n",
    "\n",
    "# lens light\n",
    "kwargs_lens_light_mag_g = [{'magnitude': 15, 'R_sersic': 1.6, 'n_sersic': 4, 'e1': 0.1, 'e2': -0.1, 'center_x': 0, 'center_y': 0}]\n",
    "# source light\n",
    "kwargs_source_mag_g = [{'magnitude': 19, 'image': ngc_data_resized, 'center_x': 0, 'center_y': 0, 'phi_G': 0, 'scale': 0.005}]\n",
    "# point source\n",
    "kwargs_ps_mag_g = [{'magnitude': 18, 'ra_source': 0.03, 'dec_source': 0}]\n",
    "\n",
    "\n",
    "\n",
    "# and now we define the colors of the other two bands\n",
    "\n",
    "# r-band\n",
    "g_r_source = 1  # color mag_g - mag_r for source\n",
    "g_r_lens = -1  # color mag_g - mag_r for lens light\n",
    "g_r_ps = 0\n",
    "kwargs_lens_light_mag_r = copy.deepcopy(kwargs_lens_light_mag_g)\n",
    "kwargs_lens_light_mag_r[0]['magnitude'] -= g_r_lens\n",
    "\n",
    "kwargs_source_mag_r = copy.deepcopy(kwargs_source_mag_g)\n",
    "kwargs_source_mag_r[0]['magnitude'] -= g_r_source\n",
    "\n",
    "kwargs_ps_mag_r = copy.deepcopy(kwargs_ps_mag_g)\n",
    "kwargs_ps_mag_r[0]['magnitude'] -= g_r_ps\n",
    "\n",
    "\n",
    "# i-band\n",
    "g_i_source = 2\n",
    "g_i_lens = -2\n",
    "g_i_ps = 0\n",
    "kwargs_lens_light_mag_i = copy.deepcopy(kwargs_lens_light_mag_g)\n",
    "kwargs_lens_light_mag_i[0]['magnitude'] -= g_i_lens\n",
    "\n",
    "kwargs_source_mag_i = copy.deepcopy(kwargs_source_mag_g)\n",
    "kwargs_source_mag_i[0]['magnitude'] -= g_i_source\n",
    "\n",
    "kwargs_ps_mag_i = copy.deepcopy(kwargs_ps_mag_g)\n",
    "kwargs_ps_mag_i[0]['magnitude'] -= g_i_ps\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## setting lens model parameters\n",
    "And finally we need a lens model. The default lensing units are in observed deflection angle (arc seconds) mapping the image to the source plane (reduced angles). In a single lens and single source plane model, this is all we need to specify and no futher cosmology is required."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "kwargs_lens = [\n",
    "    {'theta_E': 1.2, 'e1': 0.4, 'e2': -0.1, 'center_x': 0, 'center_y': 0},  # SIE model\n",
    "    {'gamma1': 0.03, 'gamma2': 0.01, 'ra_0': 0, 'dec_0': 0}  # SHEAR model\n",
    "]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## simluate image\n",
    "Finally we can simulate the images with the ImageModel class instance and the lenstronomy parameters for the different bands. Note that in the specific example we included a point source (representing a quasar) in the center of the lensed galaxy. The SimulationAPI supports various options. Further down below we simulate multi-lens multi-source plane configurations too."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA80AAAFTCAYAAAADRe7vAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAABnU0lEQVR4nO3deZRkdX3///etqu7qfZmenn1hHxhgEBlBQAQVRDTExG8iRzFBNPlGg1Hka46Z5ITlq3HgaDwE428Eo+KJ8sVIRNEEEFk1EdkcA7LDwAyz9Uzvay333t8fOE13vd+37ud2V3d1Dc+HZ85xPrzr3k9tn6l77+e+Pl4YhqEAAAAAAAAlVe0OAAAAAACwUHHQDAAAAABABA6aAQAAAACIwEEzAAAAAAAROGgGAAAAACACB80AAAAAAETgoBkAAAAAgAgcNAMAAAAAEIGDZgAAAAAAInDQDAAAAABABA6aMSduvPFG8Txv8k9DQ4OsWLFCzj33XLnuuutkeHh4Wv2VV145rb70z549eyZr9+3bJ5/61Kfk6KOPlsbGRlmyZImcfPLJ8tnPflZGRkbm+6kCQKwDY+IjjzwSWZNkbPvxj38sZ555pixZskSamprksMMOk/e///1yxx13iIjIWWedVXZMPfDnyiuvnMunDQDqN2Emk5GVK1fKhz/8Ydm5c2e1uwc4yVS7Azi4/d//+3/l0EMPlUKhIHv27JH77rtPLr30Uvnyl78st912m2zYsGFa/ZYtW6SlpUVtp6OjQ0RE+vr6ZOPGjTI0NCQf+chH5Oijj5be3l75n//5H9myZYt8/OMfNx8PAAtZkrHtS1/6kvz1X/+1nHnmmbJp0yZpamqS559/Xn72s5/JzTffLO9617vk7/7u7+TP/uzPJrf/8MMPy3XXXSd/+7d/K8ccc8xke+kYDABz5cBvwomJCXnwwQflxhtvlF/84hfyxBNPSENDQ7W7B5TFQTPm1HnnnScbN26c/PumTZvknnvukd/7vd+T3//935ennnpKGhsbJ//7H/3RH8nixYsjt/eNb3xDtm/fLv/1X/8lp5122rT/NjQ0JPX19ZV/EgAwx1zHtmKxKJ/73OfknHPOkZ/+9KdqOz09PSIics4550xrb2hokOuuu07OOeccOeuss+bmSQBAGVN/E/7Zn/2ZLF68WK655hq57bbb5P3vf3+VeweUx/RszLu3v/3t8vd///fy8ssvy3e+851Ej33hhRcknU7Lm9/8ZvXf2traOFMJoCa5jm379++XoaEhOf30083tLFmyZE77CQCVcsYZZ4jIq+PfAffcc4+cccYZ0tzcLB0dHfLe975XnnrqqWmPO3BL37PPPisf+tCHpL29Xbq7u+Xv//7vJQxD2bFjh7z3ve+VtrY2WbZsmfzjP/7jtMfn83m5/PLL5aSTTpL29nZpbm6WM844Q+69995pdS+99JJ4nidf+tKX5IYbbpDDDz9cstmsvOlNb5KHH354jl4VLFQcNKMq/uRP/kRERF0p6evrk/3790/7MzAwMPnf165dK77vy7/+67/OZ3cBYE65jm1LliyRxsZG+fGPfyx9fX3z1DsAqLyXXnpJREQ6OztFRORnP/uZnHvuudLT0yNXXnmlXHbZZfLf//3fcvrpp0/WTnXBBRdIEARy9dVXyymnnCKf//zn5dprr5VzzjlHVq5cKddcc40cccQR8pnPfEYeeOCByccNDQ3Jv/zLv8hZZ50l11xzjVx55ZWyb98+Offcc2Xr1q1qPzfddJN88YtflL/4i7+Qz3/+8/LSSy/J+973PikUCnPxsmChCoE58K1vfSsUkfDhhx+OrGlvbw9PPPHEMAzD8IorrghFxPyzbt26ycfs2bMn7O7uDkUkPProo8OPfexj4U033RQODAzM+XMCgJmKGxOTjG2XX355KCJhc3NzeN5554X/8A//ED766KNl9//9738/FJHw3nvvrcTTAQBnB8a/n/3sZ+G+ffvCHTt2hLfcckvY3d0dZrPZcMeOHWEYhuEb3vCGcMmSJWFvb+/kY3/zm9+EqVQq/NM//dPJtgO/Gf/3//7fk23FYjFctWpV6HleePXVV0+29/f3h42NjeFFF100rTaXy03rY39/f7h06dLwIx/5yGTbtm3bQhEJu7q6wr6+vsn2H/3oR6GIhD/+8Y9n/+KgZnClGVXT0tKiUrT//d//Xe66665pf771rW9N/velS5fKb37zG/nYxz4m/f398rWvfU0++MEPypIlS+Rzn/uchGE4308DAGYtydh21VVXyU033SQnnnii3HnnnfJ3f/d3ctJJJ8kb3/hGNY0RABaKs88+W7q7u2X16tXyR3/0R9Lc3Cy33XabrFq1Snbv3i1bt26VD3/4w7Jo0aLJx2zYsEHOOecc+c///E+1valhh+l0WjZu3ChhGMpHP/rRyfaOjg5Zt26dvPjii9NqD+REBEEgfX19UiwWZePGjfLYY4+p/VxwwQWTV8NFXptWPnWbOPhx0IyqGRkZkdbW1mltb33rW+Xss8+e9ufUU0+dVrN8+XLZsmWL7N69W5555hm57rrrpLu7Wy6//HL5xje+MZ9PAQAqJsnY9oEPfEB+/vOfS39/v/z0pz+VD37wg/LrX/9azj//fJmYmKjSMwCAaF/96lflrrvukltuuUXe/e53y/79+yWbzYqIyMsvvywiIuvWrVOPO+aYY2T//v0yOjo6rX3NmjXT/t7e3i4NDQ0qULa9vV36+/untX3729+WDRs2SENDg3R1dUl3d7f8x3/8hwwODqr9l+7nwAF06TZxcOOgGVXxyiuvyODgoBxxxBEz3obneXLUUUfJX/3VX8kDDzwgqVRKvvvd71awlwAw/5KMbW1tbXLOOefId7/7XbnooovkhRdekF/96lfz3GMAiHfyySfL2WefLf/rf/0vue222+S4446TD37wg2odelfpdNqpTUSmzdb5zne+Ix/+8Ifl8MMPl2984xtyxx13yF133SVvf/vbJQiCGW0TBz8OmlEVB8Juzj333Ips77DDDpPOzk7ZvXt3RbYHAAtBkrHtwFIujIMAFrp0Oi2bN2+WXbt2yT//8z/L2rVrRUTkmWeeUbVPP/20LF68WJqbmyuy71tuuUUOO+ww+cEPfiB/8id/Iueee66cffbZzNJBWRw0Y97dc8898rnPfU4OPfRQufDCCxM99le/+pWaniMi8tBDD0lvb685rQcAFjrXsW1sbEx++ctfmtu4/fbbRcSe3ggAC81ZZ50lJ598slx77bXS2dkpb3jDG+Tb3/72tFVTnnjiCfnpT38q7373uyu23wNXjqdeKf7Vr34VObYCIiKZancAB7fbb79dnn76aSkWi7J3716555575K677pK1a9fKbbfdptZVvuWWW6SlpUVt55xzzpGlS5fKv/7rv8p3v/td+cM//EM56aSTpL6+Xp566in55je/KQ0NDfK3f/u38/XUACCxb37zm3LHHXeo9m3btskPfvCD2LFtbGxMTjvtNHnzm98s73rXu2T16tUyMDAgP/zhD+XnP/+5/MEf/IGceOKJ8/20AGBG/vqv/1r++I//WG688Ub54he/KOedd56ceuqp8tGPflTGx8flK1/5irS3t8uVV15ZsX3+3u/93uR4+573vEe2bdsmX/va12T9+vUzniqOgx8HzZhTl19+uYiI1NfXy6JFi+T444+Xa6+9Vi6++GIVAiYi8vGPf9zczr333itLly6Vv/iLv5Cmpia5++675Uc/+pEMDQ1Jd3e3vPOd75RNmzbxYxHAgrZlyxaz/YEHHpCurq7Ysa2jo0O+/vWvy3/8x3/It771LdmzZ4+k02lZt26dfPGLX5RPfvKT8/l0AGBW3ve+98nhhx8uX/rSl+SZZ56RO+64Q6644gq5/PLLpa6uTs4880y55ppr5NBDD63YPj/84Q/Lnj175Prrr5c777xT1q9fL9/5znfk+9//vtx3330V2w8OLl7IXewAAAAAAJi4pxkAAAAAgAgcNAMAAAAAEIGDZgAAAAAAIiQ+aN65c6d86EMfkq6uLmlsbJTjjz9eHnnkkbnoGwDUFMZHALAxPgKoZYnSs/v7++X000+Xt73tbXL77bdLd3e3PPfcc9LZ2TlX/QOAmsD4CAA2xkcAtS5Revbf/M3fyH/913/Jz3/+87nsEwDUHMZHALAxPgKodYkOmtevXy/nnnuuvPLKK3L//ffLypUr5S//8i/lz//8zyMfk8vlJJfLTf49CALp6+uTrq4u8Txvdr0H8LoThqEMDw/LihUrJJVaOLEMjI8Aqo3xEQCizWqMDBPIZrNhNpsNN23aFD722GPh9ddfHzY0NIQ33nhj5GOuuOKKUET4wx/+8Keif3bs2JFk+JpzjI/84Q9/Fsofxkf+8Ic//In+M5MxMtGV5vr6etm4caP893//92TbJz/5SXn44Yfll7/8pfmY0jOFg4ODsmbNGlna0SApxzOFnpfsTIDnOT+lV+sTVUviM5zJ6hP2fY7PtibfetJHJHu+yapFEny8f1efbPtBwvok/Um87WTlEoZBwvpk723y/rjVBWEo/cMTMjAwIO3t7Qn3MncqOT6u7owZHx3eCs+hyHX4qNQw47Ydh35XZjMzGN9mx+kjnvSLM4vNJB0fy2/LocZpO/FV87ovl+041Djvb5ZvSRCG0jOUP6jHxxUdjdPGR2ustL7b1phoj0lWnVFo/jadxT7MMqtOt5n/XhhN1sfLq9SgU46X1k3WVT6jyfxOBJ5DkbEt87nqttD8Aeb4OqWNFz5ltJU+BxEJfb2PILB+qxl9MV4D6/mG5ifN/mS4lAXGa2XvVz8Pe0w06sztWW3G1kpe0yAMZf8Mf0MmCgJbvny5rF+/flrbMcccI//+7/8e+ZhsNivZbFa1pzxPUtaHyJD8IDVR+QI7aE6m9g+ak5mrA7GZ1iftUKLtz3Hf7YGzkvXJJK1faNPz5nV8XIAHzS6b4qC5UkWV2UylDj4ruS2n7dTgvlzrKnX4ctCPj1MPmo2xcj4Omu3X2HUf+ujw4D5oNp7vbA6aS59cpQ+aHetM1r/d5r/nxkGz8TzsQ9yFc9Bsv37Wga/1fK3NGXVOe43oXsRQOJMxMtEl3NNPP12eeeaZaW3PPvusrF27NvGOAeBgwvgIADbGRwC1LtGV5k9/+tNy2mmnyRe+8AV5//vfLw899JDccMMNcsMNN8xV/wCgJjA+AoCtkuNjWqZfac5YV5qdb/+r7DRp84KvdZXVeZ6L8VjHfVhXn+27rIyrm9aFTKvO3JrbVHYvpadsm1dkrc6UXsg0ZjDbU3/dpvRKym2KsPlmmM/B7UpzYPTZmMVtXsl1vTLseqedPXPa6p/RZn5+3K6sm58ftwvr9ntewUk3ia40v+lNb5Jbb71V/t//+39y3HHHyec+9zm59tpr5cILL6xcjwCgBjE+AoCN8RFArUsUBFYJQ0ND0t7eLss7GxPc00wQWKX6khRBYOUtqCCwxPW1GwTWOzQug4OD0tbWlnAvC9uB8XHtopjxkXuaK7EZ7mmuUBDWq9tyqHHazgILAqvofd/zEwS2ZzB3cI+Pnc3TxseMEbw0P1eaZ34VmCvNC/xKs/H7qOJXmn3jSrMVBBb6Rl8cf7/Nx5Vm47Wy2sz7q11fZ+uKtGObX9LnV8fImf2GXDiL+AEAAAAAsMAkuqcZADD3PK/8VdlKTTBxnOzjdNWmYleaK3UV3WFXFeWyQ5eriA7bqVQytMu+XGe8uFy9cLuyHb8hO9m2tCZeUKk3zfXy8Dx8KBdWZjYAHDw4aAYAAEBN8H73v8m/m8tB2Y+beZ3r8lKzWa7K2IPrklhWlXFW1F6aym3atbkUkHPwkts0Zs9a49i6RbO0z1Z2lDk72215JAn1Ps1TY8YMc7GW0rJ2Yb1OKWuqvLUEk9U/t6WarFtYzdfKmCoeWp8V6611DO4yl90yqsyp3cZ+rSWtS79rszmxyPRsAAAAAAAiVO1Kc8pL2We8LEmDvRLOXUwl3r51aqlMvesN+yLimaetytQnPO3hHj7xu/qEp2Sc39MD27fOyJV/RKJql2l8U/l+snAs3y8mq7fWDojgJUyEsYIXykka7JU0mKxSgUalmH4IAACA+cSVZgAAAAAAInDQDAAAAABABILAAAAAUBNSKW/aOs1mnpL1wFkFhjmGb80mCMyxL/aa0Y7PzXqsua6wbrITmhyvvbnlqJl9MW/jK2my7k4z26zb3lxDv6xW645K8zZF3Rnf6EvKCsay+mc+DbdQNvNWO70UtPk0zJsAjVv9rM+ZfQOh45rMjjfmmd+Dkn0QBAYAAAAAwBzgoBkAAAAAgAhMzwaABcbzvLKrALiE1LtMQXJdacAlFd9a+1HVOPTK6bm59Cd+M05FlU1rd3pXYisqFUzvEtAfOKb4u9SZUyNnUOOS5O/WbYf+VPD9SLa+wcywugAAzA2uNAMAAAAAEIErzQAAAKgNnky7pG4HbRkPcw3Gcg4Cs0KrjH24bs/xeVizH9yfr9FmBG25Pg9zRocRDGUyQ8mMvqR12pZXso/AM+ZxGLOfrJksZqiWNX/E8bWzBMY8kzAw2lJGm9U/67kZQWDW25PyjTrj9bPeRjOozNhLaDw4ZfXZMcavNMxLxHx7zc9B6ceMIDAAAAAAAOYAB80AAAAAAETgoBkAAAAAgAhVu6fZS9n3TJi15urh0VLWAt/ltu/akQPbT5hd6nrPg4hIKmlfrEXpy9Yn234mU5+ovr4+WX1dNpuoPp1O9pH1jXtGyslNTCSqnxhPVp/P551ri8VCom37CaNZXVNxZ8q8J6hcPbGvAAAAWIAIAgMAAEBN8LzUtGXwrGsHrqFa5nJ6ZgCS2RHjoVYAkmPol/FEXJcFdA05s56c/VrpZ5yyrl9Z58YD/VjzFLq137S1X32oUvo6B8Y6dKFfNPrhFpZlNYbme2a8KOYFCSv0y+ifcfXADP2yAuKsMC8r9MsMPrP261hntFnXcErD20REAvPjrTtoff+s/nnGBkuXw3RYHTMS07MBAAAAAIjAlWYAWGA8z7OvgEz+d5dtxNeU20fybcWfg3XajlN/4jfk9Bq5LD7h/Bo5PH/zUs101lWVmWzHpT++wy0aflFfDbG43E7isq0g8GNrvDD+XhSXO3Nc3tnA4TYTe9kUszC+xOEySLm3jbtcAGBucKUZAAAAAIAIHDQDAAAAABCB6dkAAACoCZ7nTb9Fwwzk0szAMCsYywzLsvsx0zbrdhYzCEzv1t6HUWc1zqrP1gtotVkJTWbik9GWNvZrtKVL9hsUrTAqfYgTWLc/WE/BKjNfE10XOIZvSdotWc1+OY3QL+OWlMAKILP2aryPodHp0Hgf7dtTrMA143NmvtBOmzMLS0O/ftcYv31HXGkGAAAAACACB80AAAAAAETgoBkAAAAAgAgcNAMAAAAAEIEgMAAAANQET6Zn+VjrzTuHfpnZVrOpc2yzQr+ct6f74h5eZgSQmaFkVp1bSJe5Rnyo20KjzQppShlHKunS/Rqvp2+GgxlBVlZ3jeSptPWCGnW+GSxm7NfMULMCtIx9GHlX1hL3Vi5WkLISw4zXyngfQyOAzHpNrc+j9VqlrPfDeEPsPusmT3S4WsoqnKGqHTSnvJQ54Ni1ybZtDUZxfUlWn2z76stdrjZlpemVqc9kE9XXNzQmqm9sbklU39q5KFF9c0ey+rqE/S/k84nqR/r7EtUP7t+XqH54cMC5dnx8NNG2JZ9LVu8ZI2w5CccdMy2yXL3rDoKEGwYAAABmgSvNALDApLzyJwtdztu5nNxzPf/nsi1zORK1ncrsy7riMaPtOJykTKXrYmtERFIOJzDrsvEn/Rqam+NrmppiazJ18f0uFIqxNaOj47E1IiJjI8OxNRNjI/F9ysXvz/fjT4Z6DifhnM6/OdQ4nyB0KLSuSKmacp3inCIAzAnuaQYAAAAAIAIHzQAAAAAARGB6NgAAAGqC53nTbtFwDeQy66ztG/eR2MFi+rqTdZuKHb5lhWq51dmhX65BYMZ+retnRp/NwDWrLxkjjMmzboXRhyDm7Qkpfc9BJjO9f16d3r5f1I8LjNs2zBs5jOeaNl6n0LzloqCbQr0X63MRhvqWmWJgbK9ofFaMnphZa0aX7bs6rFAyo8podMztitqi3p712ls7sYLPSj5TnvUCOOJKMwAAAAAAERIdNF955ZWTZ/gO/Dn66KPnqm8AUDMYHwEgGmMkgFqWeHr2scceKz/72c9e20CGGd4AIML4CADlMEYCqFWJR6tMJiPLli2bi74AQE1jfASAaIyRAGpV4oPm5557TlasWCENDQ1y6qmnyubNm2XNmjWR9blcTnK53OTfh4aGZtZTAFjgGB8BIFqSMTJqfPS86YFWVjBWyjV8yzUsaxZhY1ZfUimrf0Zb2jH0S7W4vwbm62f0zwwHM+rM52uEg6VTRhCYFfikeyLp1PTtZer1uvS+sS3ft+KojD2Yr50uKxpr3Vv9zaTdQqt83wgvC6w+G++j8WE2w7eMsLXQenLGfgPjY2GGfll11n6NFDGrLbDSxoztWU8jpYLk7EgyF4kOmk855RS58cYbZd26dbJ792656qqr5IwzzpAnnnhCWltbzcds3rxZrrrqqhl3EABqQSXHx1TKM1NYD7D+YVDbcKpxKHKsSzvs0KlP1o81tS8ribVkO2n9I0ptp74xtiaTbYmtERHJNrfF1rQt6oqt6V6+JLZm8dLu2JqGpvjnNj6ei63Zu2d/bI2IyJ5XdsXW9PXsja0ZG4zfX358OLbGl/jn5gd+bI1L0GoQuKaxxn8Bwogc2+k1Zbbj9pWuiqRjJL8fASwkXmjnpTsZGBiQtWvXype//GX56Ec/atZYZwpXr14ta7pby/4onMqxbJJxgihm+8ke4PpD84C0caYwstbhx+C0+kw2UX19Q/wPqakam91+MB7Q2rkoUX1zR7L6uoT9L+TziepH+vsS1Q/u35eofnhwwLl2fHw00bYL+fgfiVO5/GCcyjzZWUbSoSV0PPsXBKHs6Z+QwcFBaWuLP1CpltmMj0ctayl7EMpBMwfNcV73B82Fyhw0m1c4SmscD5pdylzGzXIlQRDKy33jC358FIkfI6PGx3XLOiU9ZZyos652mldyF/iVZqttoV9pNpbJSmVe31eazd9i1g8o80qz3l6+oJecCorGVWDjuQXGGGeNaaGxJJZ1hTuYRZ3rVWXXK82hMaBa765fUhcEgWzb1zejMXJWCQwdHR1y1FFHyfPPPx9Zk81mJZvVB3YpL3Q+GHY9uH6tPtlBsMuPven1ybafycT/eDugLtucaNv1TR2J6pva4n+0TdWxZGmi+u7V0VNRLZ2r1yaqb2jrSFSfGx1JVN/7ysuJ6ve8EP3Zt3gJth/2xv/AnG4wWV8cflRO5cscH2S7hvnPYo29+TSb8READnZxYyTjI4CFZFbrNI+MjMgLL7wgy5cvr1R/AOCgwPgIANEYIwHUkkRXmj/zmc/I+eefL2vXrpVdu3bJFVdcIel0Wj7wgQ/MVf8AoCYwPgJAtEqNkd7v/jf591lNna7wY42Zi2lrqq9xi4k9ZdtxerY1ddr5+TpOZTf67KWN52FMKbembNcZy42lzCnberZm6etSV6/7Eab1tgJft3mBEVCVsgK+9Gy7fFq3paxl1Hx9q6BvTLv2Q/06WdPYPSNYzDP6Z00p98ypzsZnwJrsbE0hTBn7NUK63KdY67aU0RamdZtnbK/01r+ymRAxEh00v/LKK/KBD3xAent7pbu7W97ylrfIgw8+KN3d8fdXAcDBjPERAKIxRgKoZYkOmm+++ea56gcA1DTGRwCIxhgJoJbN6p5mAAAAAAAOZhw0AwAAAAAQYVZLTgEAKi/llV/T2AppsbYRx3W5vbTL/hy25bJcX8YIllHbqYtfg7iuMX79xfrW+CX4GtoWx9aIiLR1xd+XuWLNytiao9YfEVtz+NGHxfensz22ZqA/fr3jZ596MbZGROTJ3zwVW+M9+2z8hozQG2NLsRUFGYrfTDF+2T3Pd1hqz/Xyg8MyfC4L6pXbjOPS6zXN80Lxpiw96BnLEHrGK+keeDWb8C2jv1Y4mNFmrnlstJnhYNZa1eb6y27rSJvPzQwCM8K8zLWbdV2mTrfVZ4wlaq3xvm56X+qMdaDrG43lyjL1qskL9GOD0G395fGcUVfQoV/FsTHVNhGOqrZ0oB8rovfheVab8bkwQs5SxhrKvrWuspX5ZQ0+VviWuRy2EfplfHet/lnbM4PFjN2qr8EsxkiuNAMAAAAAEIGDZgAAAAAAInDQDAAAAABABA6aAQAAAACIULUgsFTKcwqOOVCbhGu4zWS9EVpQTl2mLlF9fWOrc21j29JE227qWpOovn3x6kT1i1cn2/7Ko45KVL903fpE9Y1dbqE8B0wMDSaqb3r6iUT1gbQkqi8U40OOJmsLDqkxU/sS6GCIcsLQIeBm+h6SlScMWwgcN/96CLoBANi81Kt/pv5d1ViBXFZgmBWeZGwvZQWGGfuwQrXMYDFrH0ab9XvWy1j7cNuv1b+0EfBlbs8IaUyn9e/hlNHm1TXox9brtvoG3dbYpH9n1bc0Tft7Q7P+nd3arn8vNrZ2qjYR/bzyuXHVNjqsf0+Ojui2oeEB1TbY26vailbo5bgODPMLE0abfqhv/AYMfCN8y/qtZfwcDI2AL6vO+rIFxg9Az9ie9Z0047ys765V5rQ9l7hFG1eaAQAAAACIwEEzAAAAAAAROGgGAAAAACBC1e5pBgDYPM8z77U7wCW2wSXbIe14g3g67bAt64a8Epl0/D85dfWNsTXZlkWxNY2LVsbWNHXHZza0dsdvR0Ska+my2JrDj4zf30lvPi625tjuZNkXkRyeWotjDsZEEP++jebit1MoxGcthL5DfoNLxoNDroMfOtz/5hrIkIrflsvuZnFLHgBghjhoBgAAQE3wZHrOpHVKzzPOLNjhYNZjjRAjK/TLDAyz2nSjdVLTDhszHmvVmduzgrus0C+rTh8eZOp0wFemvtlo02FemUZdV9/Sptoa2vUJ0RYj0Kt90fRwsM4V3apm+aFHq7Yj1h+mt2+83y/t71dtLz+/Q7UN7t+t2nZt26Pa6rZvV23Dva+otokhHRiWHx1RbYWUDgcrFo3AME8nhhWK+mShdd4vFRqNoX6trDIvMD6Pjif7XCNnXTNhS+tmkyXL9GwAAAAAACJw0AwAAAAAQAQOmgEAAAAAiMBBMwAAAAAAEQgCAwAAQE3QQWCuoVp6W85tVtCWkQTmGuZlh4Pp/dptujHj6Z/zKdfQLzMwrF7vo75VtWXbunRbS7tqa+7SdW1LV6i29qU60Kt1UYdqa2lvmvb37kXLVc0ph+k2V52LO1VbIDq4rKVLh5Q1LFqi27r09vZt06/ncI8OFhvcrUPEgmDQaFNN4ns6fStlRG2lrYCviIg9pzZjv0Y2mB0OZhQG5pIBVtifUVfaZtU44kozAAAAAAARqnal2Ut55pk7i3VGrhyX9UmnqjPOxpWTbdDR+eU0teuzTlFalh6RaNvty3Skfjmdyw9PVL9s5apE9SsP03H+5Sw/0v21ERFp0CfmyhoZ02ctyymMH5uofmifXgqgnMHeYefa0aGhRNvOT7hvW0TEL+jlCcoK4tc0ncZah6Acz7F+NusFAAAAAAkxPRsAFpiUZ08RfO2/x585cKpxPMHociIy43Dysd5Yv7NUQ4uexlaquXt1bE3ryqNia9pXrYuvWbomtkZEpKtbT9MrtfKwpbE1h3TH18ynNcvcTmw+vzr+hGnnrr7YmqHB+JOF+fH4k5VBfiy2Rny9hqnidPKvcguQJjzVqBizLAEAFcD0bAAAAAAAInClGQAAADXB+93/Xvu7VePYZmUYWcFdZtiYfqw1K8ea0ZM2g8D0dSwz9MtqyxhhXhn9Ez9tBIFlUnW6zbgXrq5V3+7WvGSlautYrWfnLFqtZwctWrVM9y+jA8gasnqG0qIl02fArO9uVDWVduRi/drtWKzDzKShRTWNF/X7nff1LW917fr2z1yhqB+b17Nk/KJuSxttoRhtRjhW2vFLFBqXX83JaWbAl27zzIAvt32EVg5YzN+T4EozAAAAAAAROGgGAAAAACACB80AAAAAAETgoBkAAAAAgAgEgQEAAKAmWcFdnhUwZAV8WQFDxj6sNnO/VuiXGfBlBHc5BoGl0/qne8YI/arL6ICvdJ0O2qqr18FVDW166bvWlYeotsVH6aX9Vh27Xtet1oFh6VRWtWXq9evS0qpfg7bSGlVReaX7FBFpN9rqV+vKhq43qLbQ16lVOwb1cnotbXovYa5L79hYQm8i0CFiYWAEgYU6lCw0ltsLQ/3+hNaSe1YilxXB5eu2cI6/z44LBJq40gwAAAAAQASuNAPAAuN5nnkW9YBUmf/2Wk38fqylUCzW1Y5S1lWNUtnG+OsBzZ16GZJSHSsOj61ZdIi+2lGqc83RsTWtXfH9ERFpb4lf8qSxQV9ZKbXQzmS7npWvq3d4/q2LYmua2xfH1oy39sTWFEb6Y2uC3FhsTRjkY2vEuCJjlolDndOmot+VcDbrqQAAIlXtoLl0nb1yUg4/2Kaypr2UU1cX/0NmqoaWjkT1bd16zboonauPSbTtrlXHJ6rv7tZr5ZWzYrE1ISXaaocfjlMtSfgJTPbOJv+A9zv88J+qqbEjWX2rMaUmQmPCz9nEoF7fr5zixGCi+tB3+PE4td6cnlOGw4FggjIAAACgIhbaSW0AAAAAABYMpmcDAACgJpTevmJNPrLDwYw6x4Ahz9PXmFLGdae0VZd2C/jKGAFf6ZRj6FedniVXn23QdVk9I62hTd+C0rZS3wKz/PhjVduhb3yjalt1SLdqs+Yg6ugpkXFjglrRuIui9OaLuiZdo3sxO88abUPWHR7GkVXWeAHWnqBD1AZ271Ztoz27VFtDs77Vye/QM0MDIxwstMLBQt1mh3kZ9494xhzQwAj4MvZhzRw0r+Y6zjAMPd2/0nA+11nOFq40AwAAAAAQgYNmAAAAAAAizOqg+eqrrxbP8+TSSy+tUHcA4ODA+AgANsZHALVmxgfNDz/8sFx//fWyYcOGSvYHAGoe4yMA2BgfAdSiGQWBjYyMyIUXXihf//rX5fOf/3yl+wQANYvxEQBslRgfVRCYnQRmNOnrRFboV8q5zQoH06FIVl06o+vSRhBYJl2v2uqMpTHr6vXSqdlGI/Srdalqa12+VrUtO2adajv8pI2q7ZC1ehlNK/QrZ7QNGyFauYEh1VYcHVFt+eHpdWPNOgRrbIUOOFvbqt8LK5DsmUHd2rtPdzhXGFZtRSNoqq5NB3cFKR1a1dii37N0Rn8G0p7+rLR0LFJt+dFR1eZP6OcRpvVnylw21PgehIGuC4vGY0MrMMyoMy7nml2x2qxgv5i/JzGjg+ZLLrlE3vOe98jZZ58dO+jlcjnJ5V77ugwN6S8EABwsKjE+elJ+YHdZq9r6MVgq7bjodToVv0J6pl4ntZZqaOmMrWntXhlbs3iVTnYttWTtkbE1Hcvj161vbHBbe74hGI+tCfsHYmt6d5XmwmqrjB+DWny/jRxUZffOvQ5VImOD8f+2pxx+rmQbjRjcEg1NrbE1OeOAoVRxTB9olPIL8T+TAt9InjWkXdaut34J6qLI/xIu8HXs+f0IoFYlPmi++eab5bHHHpOHH37YqX7z5s1y1VVXJe4YANQaxkcAsDE+Aqhlie5p3rFjh3zqU5+S7373u9LQEH9VQURk06ZNMjg4OPlnx44dM+ooACxkjI8AYGN8BFDrEl1pfvTRR6Wnp0feOGUxc9/35YEHHpB//ud/llwuJ+n09Gl82WxWstn4KVAAUMsYHwHAxvgIoNYlOmh+xzveIY8//vi0tosvvliOPvpo+exnP6sGvHJKgxzK1ybppUgqnewBdVm3s54HNLXp4INy2pYd4lzbvXZ9om0vX3N0ovrVHTosoZxDmycS1R/SsD9R/aKCDikoJx+63V94QKY3n6i+dzT+fsKpWlIu95+9ptnhfr0Dmpp1cEQ5443JXpvCqA6XKMcvxN+zOVUQJPziBo4TX5zu+Zt/lRwfAeBgUtnxsST1wdOPtX5fmiFBZjiY0ZYygsCsNiv0y+ifVZcxXoO6OqtN/3SvN4PAOlSblRmx9KijVNvhG9+k2o40Qr8sOrZLZKxPR4H1v/C8apvYuVO1+cODum18ett4Wr+eA2n9e+ulNv0bOMjoY4Ag1EFgxUb9myls1tkJE74O1Srm9G/ddIveb+eyDtXWY4SD5Qb0PpoW6+0N79e/yXMpI8DO+jwadUUjCSNtfK9C47sRGr8Jre+QFcxgbs+YLJ0yciM8CUv+PnOJDppbW1vluOOOm9bW3NwsXV1dqh0AXk8YHwHAxvgIoNbNeJ1mAAAAAAAOdjNacmqq++67rwLdAICDD+MjANgYHwHUEq40AwAAAAAQYdZXmgEAleV5MQGIDkkWLgGKriGLKYeQnrr6+CC6prZFsTWLlq2KrVl+yGGxNSsOWRNb097ZGVuTyY3F1oiIeH09sTV1w9tja0b6rBid6XbvjH+tvUYdHlOqbyL+J8DuPrfgvdxgfF19WIitaaqP79N4Y3x454RDwGe+TofplPIdPvuBFWZjCI2AG1Xj8uUOol9rh13UvJTnSWrK4GUGfBnXhLyUERhm1RmhSFZQUsoK+HIMWTLbjM9Rxgi2tYLA6ur1972xVQd3da5ardpWH7dBtR162BLVZn0Txoq6bXS3Hgt7n3hatY2//IxqS43s022+EaIVlASTGgFQQUq/ThM9+nXPh7oun9KhX0GTEaTb1KGa/Gb970o6p9+LbKjrjNw3qbM+U4EOuR3bo8PW0p4O7soYr0tofB4t1tAT6F1I2hyIjDZjv6ER9mr1LjRCyaztpUoCe2cTJcuVZgAAAAAAInDQDAAAAABABA6aAQAAAACIwEEzAAAAAAARCAIDAABAbfBS0wJ/PCPR0DPCk+w63WYFcs1me1ZdyojV8ow2K2ysLpNVbVkjpKpl8TLVtvTwI1Tb2mN1sKJ1cDBh5C7ljNCv0cefUG3BtsdVW+PIbtXWEOrQr0zaV21hOD2BzAqPCtJG6FfBSNoq6GcbpPT2coUB1eYP9au2YosOcyzmJ/Q+irqtaH1WfB2gaAVtFYx91GeMkLO0EQRmfG6tUMLA2K8VcGjkspkRh1abGfplFFr5i1afS8P+ZpOVyJVmAAAAAAAiVO1Kc+ySKlOkXAsP1BvLCpSTqYtfmmKqhtb4ZVOmaluiI/6jLF59eKJtr1p7SKL6o9qTvZbHt+1JVC/d4/E10+iziuXUWxHzZRzWb6yHUMZQJtlnp7ctfumXqQYX62UHoowNuteKiIwPtCeqzw3vT1RfyCV7r/xistfeOlNseh0sqQIAAICFgyvNAAAAAABE4KAZAAAAAIAIBIEBQI1xmaHuOVRZATUWl1te6urjb3NpauuIrelapsNrSq06JP6Wl0MPXxVb097k8Pz73G5LKA72xtZkci/F1mSH98bWeAPxtzKMSfytI2NhZ2yNpLrja0Skra4jtiZsd7gVqtAaW1IcaomtmWiM31euri6+O0YAVKm0lUhjCBzuLgodvpNWKM4BCe9mO3jNJnXIqPOMOs8Ii7LqrM+H9bGybkVMGWFWGWOsbWzpUG1ty/VYuvxIfQvgYmN4t26uKo7pNr9HB4Flererto5in2qrr9dja1NK3+LXkDbGOytpqkQh1AFiYykdqjVkvO7jxnd11AgRG/V1P3wj0C0wPgMF45a0ghH8lk3rztRljO2pFpG6Ov35qTOO/Hyjf0Xrc+sYyGWHdOnna72LdpiXG/M2v9IHz2KM5EozAAAAAAAROGgGAAAAACACB80AAAAAAETgoBkAAAAAgAgEgQEAAKAmhGEo4ZQgKDNgyEhdMwPUjEZzew7BUyIinhVtZIV+GeFbabNNN6bqdMhftrlNtXWtWq7aVh+yVNfp3ZqhUhNG/l5gBJWl0/rRjVkdLdaa1q99S51+vvV1xmtaGoRlvbkF/biRcd2PvgkdGDZohH6lc/q5BqI/Z56f03Whfk08X/clndWHZX5bk2obbdJhcPkxI4DMSDRLGyFinvFepIw0L89sU01iZOSZQYieFdyVMgqN75/1WJc2c5+OuNIMAAAAAEAEDpoBAAAAAIjAQTMAAAAAABGqd09zGNqrWleAl3Dl6lQ62ctQl9X3EpTT0NzqXNva0ZFo24u7kj3XtYsTlYs0GzexlJX0I5X0vI1e+L6slflE5WtC9/dKRGRfQ7LnO5jW96ZEGR4fSLTt4d6diepH+3Ynqs+NJnuuReNep3Jc7xkz70uDweX1dHsxPeumpRKpTF1sTbZR34tXqq1D35tXasnSRbE1qxzGxs7YCpGgO+tQJeL3G/dhlagfm4jfUG7EYWf6XrhSo378+5+u646taVq0JL4/IlLoXBlbMxLGjyFtu+LH4LAwGlsz1h8/vo321cfW5BzGsaLn9j1KOdQFDvfbGbcVTmJ8BIC5QRAYAAAAakTJCSo74ctock4n0m1GspF1otfaheP5YPE8I/QrrU9GZowLPfWN+uReS2e7anM9bW+dArXaWg/TMWJ79nTo/Qb6YlNbWgdmSYvxfjQbCWlNJSdgM8YJsBF9krKlb0y3Det+9IzrZ5uq1yd9vUC3TWT1axJ26pOPYUuHaivW6fexoU1fvBob1CdX82O6LSjq1yDToF+r9Ljx2UsZIW/GOcTAaguMQC7jbbSC81zbrHP+5jnD0r4QBAYAAAAAQOVx0AwAAAAAQAQOmgEAAAAAiMBBMwAAAAAAEQgCAwAAQG0IpSTnyzV9S4dMhWZ0kFFnhY0ZbWFoPNYIG7Pyx4zsJAmNa1uhcbkrVa9/zqdSuq3ii9a0LVVN3WvWqrZ0sE8/NjBWQ9EZWiLLdKCZHFqasm/FlBmrFWwd0m07B1XTkkEdlpXy9PIzzW3LVVu+Q7cVm/VKBWO+Dv0q6N7JhOhVYMaGdKDZ+IB+Hvlx/RqkMzowLJ3Wn+W0EdwV+Ead8Xm0PqPW5zuwwrwcA76sPC8zMKy0aRZfAq40AwAAAAAQgYNmAAAAAAAiMD0bABYYNfuw9L87TC9ym4HkOk/JpS6+xpw6VSLt+bE1dan4GmPVTqM/8dIuGxKRtJ5ppzU71GSNOXGlrLlvpbsKSqcvGrsyphaWKqw9Mr4/IpJavDq2ZtRhO02LFsXWjI0aUyxL9O56Obamf5deY1Xty5qjWCKVcvkk2dMO1bYcisxliRPsAwCQHFeaAQAAAACIwJVmAAAA1ITwd/977e92jWozL9EbdeZejdCviN5pRuqXy3QhkYjEMD3TJigWVdv4oBEW5RvhW+n4GRdJpNcfpRsHjJkfuV7d1m5M2Tm0w9hLadjYqUbNL3XTuueNfhifi3SDalq86HDdduyJ+rHNh6gm37hGac1hGTU+PzsG9Xvm1elH9+/aq7c30K/a6o3pU3ljRo1vzKAJrDbPuP5qTHkxw7yMVs/4npqhX8Z3yN5HWPbvSVTtoDkULyK10Kg1xozy205Yn3AHoa8Hp7L1RWOQipLXqXblpK2ovTIa53zqVvyUwOlc5jROpVMEK6lzheNczN9p1+NqWW1Bk3Nty75dibbd0NqZqL4um6zzKYdpitPqE84T9BzrmX4IAACA+cT0bAAAAAAAInDQDAAAAABAhEQHzVu2bJENGzZIW1ubtLW1yamnniq33377XPUNAGoG4yMARGOMBFDLEt3TvGrVKrn66qvlyCOPlDAM5dvf/ra8973vlV//+tdy7LHHzlUfAWDBY3wEgGiVGyM9mRb5YwbZGOEXRiCGa7aOnTDktg/r8pT9UKN/Rl3g6yCw/MSwahsd0EFbfXv1cm2rVlQ2CEzEyEtJ1+m2lJFpk7IyhqxDlTeX/L3bqDldNzW+ottarKA24zVpMp5Xc5exX7fQL0uz8UE7ul1n4eRWrlFtVq5Npk6/7l7aeD09q4fGh8/4QFoBe1ZwlxV+Z4Z5WaFfxuZMM8/4cpLooPn888+f9vd/+Id/kC1btsiDDz7Ij0IAr2uMjwAQjTESQC2bcXq27/vy/e9/X0ZHR+XUU62o91flcjnJ5XKTfx8a0me5AOBgMuvxMZSyZ0xD8zRuSY25vMp0gfOqJ/ErDPiF+GT73Fj86gAj/cYyJCX69+6Mrdm/KD7JP2yOT+/P9OllPCzZ4fh/21L5XGyNuTyN2pDDdYv6+JUAMq0t8TWLF8fvy5HLWgmjy+NXAWjriu9TY7PDc6uLf41Sqfi72NyT/13q4r+U5TZTK6sLuIyR/H4EsJAkDgJ7/PHHpaWlRbLZrHzsYx+TW2+9VdavXx9Zv3nzZmlvb5/8s3r16ll1GAAWKsZHAIiWZIxkfASwkCQ+aF63bp1s3bpVfvWrX8nHP/5xueiii+TJJ5+MrN+0aZMMDg5O/tmxY8esOgwACxXjIwBESzJGMj4CWEgST8+ur6+XI444QkRETjrpJHn44Yfln/7pn+T6668367PZrGSzLpOyAKC2MT4CQLQkY6Tr+Bg6BnyZoV9W1pG1DyOwyGgyb4sJAt3mh/o2jGKgA76KxmPzBb3fiZEx1Tba36Pa9u/ap9qeLOpbaxYvMgKk8vo6W9Y32kZ3q7a6oX7VlsqNqzZpsG5PaTbaHir5+2lGzS91U0G/xpI3+jFi3G6UG9Btu/RrLCvajb5UVmOHvv0o02TckmLcyuMZ10s9z7iGat4pYgR3eQ63FImI5xjO55zN53orSmndLG5hmfU6zUEQTLvnBADwKsZHAIjGGAmgViS60rxp0yY577zzZM2aNTI8PCw33XST3HfffXLnnXfOVf8AoCYwPgJANMZIALUs0UFzT0+P/Omf/qns3r1b2tvbZcOGDXLnnXfKOeeck3jHYRiaU1uiapMIjCku5RQdUl+nyo8lS3Ac73dLXxURGd6zLdG2BxuTTRYYSDgjv/uIZK+NuT5fWfEJp7MTn6A71UTCaRuBsfRgOZ5Dou1kbSbZtF0v7boS4O/qE8asJq1PaqYzbRaKSo6PAHCwYYwEUMsSHUF94xvfmKt+AEBNY3wEgGiVGiND8UruTzYurJjL8ul7L8NQn2y2rtMExvp8gXEvp+/rNmvZMqNMrJX9AqvPKeN+6IK+0Tk3PKDaBnZuV231qZWqbWJE30fcYlwkyIzrfeSfeUI/dugl1bYsrS9AZawX4dfGxaHjJ0oe+IyuGS3qtpeNC0EDg7ptUB8eTXj7VVuuUd/T3L7iCL29Cp/uzxjXSLLN+l7qTJ1+Hp7x3bAuilifW2sFPvO2ZKPQMz63YmQAWH0x77k2mM+t5LU37612NOt7mgEAAAAAOFglTs8GAMytuNtXrIRWaxuVqBER8f34W14K+YnYmtGhgdia/bt2xdZsf/b52BopxCd67m+OvwWiaUynzVpajKTaUl15nXBbqjljXB0p5XAnRlgf/57lWowY3hJ1ObdbdNLZxvg+OWzHuGCm+MX45xZYl/JUh5Ld+hWlkreWON0FU5luAwAS4EozAAAAAAAROGgGAAAAACAC07MBAABQE0KZPkPdul0ltDKHjAny1kz30AoRs5qM6f2hURiYdbovgdU/Y1WMTFOTavN9fV/D6KAO2urfq29/ydTp/hVa9cojo3m9nnZx107Vlup5TrXlRIdtefX61o/lxmuVyhu3rPSVrNTSYBzOFIzbOMaNthH92uVH9K07AzKs2nKNenWczJB+3ZvbdEiXK+tGk34ju6yxW98i49XpFWqKef26h9YtWObtK24hXXbYmFWnn525UosZDmbdy2JdCy7tC0FgAAAAAABUHAfNAAAAAABE4KAZAAAAAIAIHDQDAAAAABCBIDAAAADUhtLsIJfsH7Fzjcyso8AI5DKyg8wAMqMtMNqsZcR9I70sVd+s2jqOWKnaRrcPqLax4RHVVj+o69L1daotPzah+zKqt+fv00Fg9SO6Lkjr8Ckv0M/XD/Qb0jqh+9I+PL0tlbEWrzc+GL5+L0byuh8D+XrVNpwyPixGONrI2Jhq2z+i99G0olO1deg9GBFqIjt/+5Rq2/GIbsuP6fAy3whW84v6/QmM9ycw3h8zdM8tS0/sUC7H0C8zLyz+sXaAmJuqHTSHYWgPVgYrebAc309WXyjoD0s54yPWRzja0L7dzrW9LS8m2vZuI3munG15nbpYTnfSyQhH6MTB8nTaX3kJtz+U7MuxY6wjUf1+/W9DWUP79zvXjg3sS7Tt/KhObCzHL+jBvpwwMNIVyz8iWbmVWDqbuhoWhDFP0+ElMP5tU6wfKHZd/DhTyOsE0lJjxo+pUr374r8jdS/rH2ql8gXrR9R0nS3x42FbOBpbIyLSWYz/fuQzOslU7U/0D9hSqWL8uDxebIityQ/E9zmzeyC2RkSkfkV8Ouy4w3C886We2JrePTqxttToUPy/04Wc/jFeymXcsxKTAQAHF6ZnAwAAAAAQgYNmAAAAAAAicNAMAAAAAEAEgsAAAABQm6xbys3kLrfQocC4nuSFRptRZ0XqWNEzvnGrvF/UfSmM6Xvze57eodqsjKqgUWc25CZ0hk8up/NNQiPvIusbuRWezo4ohjqXIVXQ4VOplG4bNzIEBou6L/tLAhLq0romk9EvvLEpmcjr55D39OFRPqvrMmn9XIdGdBDY8KjOTxir15k+/mKdRbFnQL/u+59/Xu9398uqbWKwV7UVc+OqLTDeWysIzAr9sqJRrM+PXaib7O+ubvLMELG5zZfgSjMAAAAAABE4aAYAAAAAIAIHzQAAAAAAROCgGQAAAACACASBAcACE4ahhGXyLMr9twOszA1V41IkIr5DXcFKtikxkdchNKVGjMCUUv2DOmilVKpFh9uUCjItDtvpjK0REWnvao+tSXdYwSUl+8vGv9a5CSOQp8TgSPw58SE/vs/FffHvmYhIMLYvtmbMCBwqtXfHi7E1u17aFlsztH9vbE1+fCS2xrfCj0qELl9IEQldQmrmNsfmoBCG08dAa3iK/6b9rs4INvLMHCLdaL3toRFiZH0+AtE7toKXJsaN719fj2qqr29Wban6etVW8HX41tiY/ox77TqQKtveptrqsjoIK0jrwKzx0T5dF+hxvDHQ/bOGu7qS1zmV1q9dXUG/7kUxXhMva+xA/9uQaluq2tKLlqu24XH9b+G48SHNGJ+B3jH9fj/34COqrW+7DgLLj/Srtonh/aqtWNCvu2989gLjc2uOdcZjrXHM/L4Yr0toBpBZdY79K22axRjLlWYAAAAAACJw0AwAAAAAQISqTc8OQnG/RO44hfCAlG9MFSijYKwfV874WPz0wamGB4eca/t69fSKclpauhLVb2+Ln444bftjixPVH9avp++U09Cgp/aUMzqc7LPwcq/rJK1XbRuJnz441SsD7u+tiEjPjqecawd26uk35Yz270lUnx8fTlTvG9O6ynGdsvhafWXrAAAAgErgSjMAAAAAABEIAgMAAEBNCMQXb8pURWtuoWfNSDJCujyj0DOmQXrGXjxre8ZjfWPCm7VfX/RsrqIRRJcvGMFdGf1YK9wpb7xY2bQ+FEg161mDqbZW/dhFxpNr0QGDhSEjkGpUz9Qrig7CSns6WKtYnP5EzPfRyqdK6dmNYZN+XunWbtXWtGqN7kernu2ZD/XrWb9Yv55Fo4MvP/a0atv9zG9V22i/Dl7M5/Xr6Rd16FcQ6s9KEOrXOAiMNqvOCrozw7yMttAKIJt5W2j2Lyj79yS40gwAAAAAQAQOmgEAAAAAiMBBMwAAAAAAEThoBgAAAAAgAkFgALDABGFEkM0BDstuWSEopfyyO3lNymHZv6LDUn+FYvx2Jvz4mvEgHVuTSzXG1gQOS/Y1rXJbdm/FEfHbOrSr02FLLsvkjcZW5Hb0xdb07NGhO6X6xh26IyKj/Trop9Rw/97Ymv3bn4ut6dn+Yvy+euP3lR8fia3xjSCmUlbAjVnn8HVz+UaW287rYUm+MAynLWlohwkZIV1mnbGDwEzuMtqMBxsbNPdrfGaKRvCSZyz16BV1W11KXwMLUnqcDDzdlmnRQVgtK1aotsbmZv3YlH5d0h16zGzML9eP9XVIVX1Rb68+yKq2VEnYVnF8UNXkC3pZzXTWCD3r0qFfqdYO1RZkmnRbSvettVs/tqGjQbXtfO4V1bbrGb0s6fCeXaqtkNfPLTehx7NCXi+jGvhW6Jf1GZ153WwCvkJrFLTCu8wgMIfHEgQGAAAAAEDlcdAMAAAAAEAEDpoBAAAAAIjAQTMAAAAAABESBYFt3rxZfvCDH8jTTz8tjY2Nctppp8k111wj69atS7zjwCFY5gAv4aG97xjKcUChqG9iLydX0CEM5YyNxweJHDAykWzbQ36yLLe+rA48KGdPe7L3NuWSczNFXbJyGdS5EWXtHJ9IVL99T3yYzVR7dm1LVL/v5ceda/teeTbRtkf79iSqz4/HhwlN5RfdP8ciIkHCRBrX8oUadFPJ8READjaVGiNDCSWU137nhUbol3s4mP4HxXpsGBpBW8ZvTc8IxvJ9K9zPLdjIM0KRUsbPPivgK6zTIVWS1W1egw74almifyvWGUFYE0M6kMrInpIGTwdhtS1ZpNpWrzxMta1omVlm8fN7dDjYxPA+3ZbXr10x26Lamhbr/rZ26MfWG33pN97ul594XrXt3/6Sahsb6FVtE2M66HF80AhDm9BvRtEIkguKbmF1zuFgoRUOZoV0uYV+WeFgZmCYQ539ODeJDkfvv/9+ueSSS+TBBx+Uu+66SwqFgrzzne+U0dFkP74B4GDD+AgA0RgjAdSyRKdv7rjjjml/v/HGG2XJkiXy6KOPylvf+taKdgwAagnjIwBEY4wEUMtmtU7z4O+mAixapKcsHJDL5SSXe216wNDQ0Gx2CQA1gfERAKLFjZGMjwAWkhkfNAdBIJdeeqmcfvrpctxxx0XWbd68Wa666qqZ7gYAas5sx8cgDMUrd9uNS2yDy803jvEPnh9fmHLIkig4bKdg3FtVKu/Q72IqPjEh1dQWW9O+YnX8zkRkVZe+Z2/u6HsQS3Wsjq8JhvbG1gz2uuU89PbEZyoMOGRA9DtkOQzs3h5bMzYY3+98Lj4kw3fIPPEdgxZcch5caspVLNDIB8VljOT3I4CFZMYHzZdccok88cQT8otf/KJs3aZNm+Syyy6b/PvQ0JCsXu32IwQAahHjIwBEcxkjo8bHMJh+ciBMWWFeentWEJF5ksI682CdFDROTHqBETYmRhCYFUpmnMX0PH3SpujrDvrWPowgsCCj2wrGUxvP67bO1UtVW9o48egP63Cwpnr9Yi1euli1rWjU+52p7mXtqm1Po+5vSvRrnG3Xh0dW16zztwNG22/u/rVq2/lbHQw70quDygoTeoZFbtQKYNPBt4W8fiOLvg519QMjHMwM0zO+Q8arYIU9u27P9bvr2qa+z7M4szijg+ZPfOIT8pOf/EQeeOABWbVqVdnabDYrWSOtDwAORoyPABDNdYxkfASwkCQ6aA7DUP7qr/5Kbr31Vrnvvvvk0EMPnat+AUBNYXwEgGiMkQBqWaKD5ksuuURuuukm+dGPfiStra2yZ8+r9zC1t7dLY2MF51QAQI1hfASAaIyRAGpZonWat2zZIoODg3LWWWfJ8uXLJ/9873vfm6v+AUBNYHwEgGiMkQBqWeLp2QAAjfERAKJVaowMw+lhQWb2j5W9ZdQZeUVm6Jfn6Q1aYUd2ypC5E9XiW6sLpIy+FHWQUz6n26xf+EUj9Wt8fFy1DfePqLZsk95ea5OeIeAv122t+qFyiNFWSToGTGSkXb+Po8YLNWKE5g8Zn6nChH5vn7z/UdW2/dEH9X537VBtuSEdBJYf69P7zY2qtmJBfwaKvg74KhZ1mx/qJ+wb3wNrpQAzYM/4zLsGd7mHg+m2uV4+INGVZgAAAAAAXk9mvOTUbAVhYJ8KNIsd637HOiNYtt5h7dCp8kaEezkTOX0WL8rYuD57VM6oETFfznCypyoD6WT18SujTpdw89Kf7KWXXYP9ier37H4xUf3el59KVN+3/Wnn2qG9Lyfa9thQb6L6fN79cylin3UsJ2G5fdZwFnUAAABAJVTtoBkAYAtCEa/syQGHMwcJT1qU43nxGysWjflsJQqF+LNeuYn4kzkTY3rq4ExqCoVcbI3rud2kJwznmsvbPz4yEFvTt1dPIbTs2/FSbM3Arhdia4YdTlyO9u2OrZkYG4ytKTq8/34Q/7m2p+lqLif83Gqii7hNBADmBtOzAQAAAACIwJVmAAAA1IRQREJ5bQpIYEwHsWbqmNfgzRAxI/TLeKi1D2vWgfMNgylzg3p7RltuQt/aN9CzX2/OmBOTNhK+xod1+NTul/Ttbt2rO1VbyrgclzNehDHdZAaGzVSP0fbyHt02PrJXtxnhaL7ooK3dT+rb6HY/+RvVNtqzS7UVxvXrmRsbUG35CT1ryi/qWVtm6JfvGPDl689e0fgsW59v67thRt85BvGZoV/WFq2wP6OttH+zmYTHlWYAAAAAACJw0AwAAAAAQAQOmgEAAAAAiMBBMwAAAAAAEQgCAwAAQE0IZXrQUGgkAtmrgOm6lNFmbc8KJ7L2YW3PXLTMTAezQsSM0KZABz4VchOqbXT/Pr2Hot5e6OvO1DU0qLZ0tlH3JVyr2rxAH1pk63UA2Wh3m2o7We92xra+qJeU693+lGrLG6FnfXt7VdvYoK4b3q2Xv8sN6Ne9YCyBl58Y1nVGoFuxqAPIfCP0yw+tNv2Z8o0PrtVmhn5ZjzX2ERhtVpiXHfpltM0iRKySuNIMAAAAAEAErjQDwAIThmHElZIKc1x7wfccCo2z4aW8CWuRkenSI3oZjlKjfca6ISUGW/VVjFK9bXqplVJ7uvXVFcvLK4+JrVnrvPbM7G3f9mxsza4Xfhtb0/PiC077279ze2zNcE98zXi/vnJTKj+ir/iomlz8Z826glPKWp5F17h9Wa2rLzOrmdl/AwDMHFeaAQAAAACIULUrzUEYRqw0b0h45tTpqshUxiLh5eSNew7KmXA4K37A6ED8WfaphvZ1JKrv3RV/9WWqps5kN5gU/JWJ6r3QvNsn0uDOHYnq97wYf7Vlqh6HqzNT9W5/LlH94F63qzYiImP98VfTpspNjCSqLxb1/S/l+EGy9yrpFQ+XKywir16FBQAAAOYL07MBAABQE8JweuCPGSbkWQFfRp3ZZgUbGfdWWNdn0kaYV6gLPWt7VmiT8Vjx9YUeL2c8X994Hnl98js/prc32j+k2vp26RP5u55crtqyWX3bS1NjVrW1d7artlc6l6i2NYc3qzYpueDy8pN7VcmebbtU2+C+l1Tb2KC+lWN8TD//fF63+cYtR7lhHfrlGxfbfF8HlflW6Feg2wLjglNg3EriG3XWZ8pss4LArEAus85o0w81v3/W9szv+AxDxGZz3YXp2QAAAAAAROCgGQAAAACACBw0AwAAAAAQgYNmAAAAAAAiEAQGAACAmhCKJ6F4U/5u1JjhYFbAkNFmJRaljGAja6UWoyllBYsZfXGOKrOCkszQL71CRpDXoV+F3Lhqy40Mq7axvh7VNvCyXhGksVEHdzU36sON5iYdDlaX1iFizzXVqbZiSb5VMKLDt8aMMK/cuFE3oV+ngvE6+YF+nXxj9R1/YsKpLvCN98cI8wqNlUsC40MaGB++wPpcGG3WevRmmJdVZ3zmnYO7zFAuHWpnBfG5fndLq2az/goHzQCwwAShlB3ZXdIfrR+IqsZ1rpHxj6ymEz5L5XMOy/UN9cbXpOL/6fIc/mlM+/rHTamsPxDfHxFJjeqk1lL7D9XJsKUaG+LflN49Oim21LNPxC+H98KTL8bW7N7+SmyNiMhwr/5BXWpiaH9sTcH4UVuqmNc/XlWNH/959P345SmtBNlSrsvluW0rfjvllt1jST4AmBtMzwYAAAAAIAIHzQAAAAAAROCgGQAAAACACFW7pzlwuXHnAH3/d3lO999Ne0Ci6nx+LFG9N7zPuTadTnYeozelgwTKyYgOdyjHH92ZqL73xUWJ6j0zcSPa8L74++am6t2+PVH9/p1u9+8dMLgv/j7GqcYG3T8LVmBFOcV8LlG97yf73Dvc/jeNFVRRvt7te+taBwA4CJVmPlj/JDj+brTuAbf+5fKMtCPPCCcy92vVWX1O6TormsI3QslC4x/oMGUESAVG+FRohFQF+veEX9CZFP74gN5vQ6PeXr1qkrGU/r1bl9XhYNbvxKA4/bXyUkbAWahzDcJx/bys9IN8TrcWi/qxgfE7Kiwar3FpcpmIBNZ7ZoV+GcdLQWjs1/hQmZ9v13Aw43UPzb64hYjZj7XaVJNYX6zQ+F6ZvzormATGlWYAAAAAACJw0AwAAAAAQAQOmgEAAAAAiMBBMwAAAAAAEaoWBAYAsIWBFekxhUPITeBW5CYVX2jkkuia/ER8TTgQX2OEpZTyC+OxNcWx/tiaXL9bOGD/9qdia55b1BJbU+fwr/LoUHxIYM+u+NDEfXvia4b6B+I7JCK5sfiQSd8hRNM3wnZUjR8fgOkSNuoSKuhXaDuudUbelGKF+7z2+NdHUGJcDpgZgGS8tp6RtOV5RiCXMZ6mjLrACOmS0Lg+ZXXaKcVIxDPe49AzwpjSui1l7CQ0wnCtoCkJdDhWIHWqreAZIVo6p8vKPZNsQW/PL+pCP5j+mtan9Q6KKd0Pf8wIdDPeHt8ITCvm9T78ohWWpR8bGnWB8YG0wrLsz7IVYGdszzWQywglcw/4cqwzQ79cn69qMoPPrKbSpzubIZIrzQAAAAAAROCgGQAAAACACBw0AwAAAAAQgYNmAAAAAAAiEAQGAACAmhCGoYRTQorM4CAjuMsK8zIDvqygICv0y7ju5AVGupURNmay8sKMPC4ruCo0dptyDEqyXr+08XzTRnKX7+vHjk8Y4U7G65zJGK+pkRhWMPZR+jIXjTAzKwszXzCCwIxkKM8IUQt8vQ+/aLRZQYVG9+wALV1o5f5ZbVZwlxWWZQWQWRu06sxwMDO8zDUczOqKFRjmWmcVlvTFMbjRkvhK8wMPPCDnn3++rFixQjzPkx/+8Icz3jkAHEwYHwHAxvgIoJYlvtI8OjoqJ5xwgnzkIx+R973vfTPecRh6Elpn5AyOZa9t2/Ws3oF6iV++Ylp9PuFZiiRnNXwdkV++L/HLfExVHNmbqH5497OJ6htbmhLVqzNAMXIjyZ7vcH/8kjJTjQzGL+Uy1fhosv4UcvHL4BxQtNZoKCNwWIZnKt9lbZOp27dOZZerT7Z5569J0u3Op0qNjwBwsGF8BFDLEh80n3feeXLeeefNRV8AoKYxPgKAjfERQC2b83uac7mc5HK5yb8PDQ3N9S4BoCZEjY9BGJr3nh1g3Yena2Z+304plwkhoRdfFEr8TBrrnilVE8TPDirm42d1FEYHYmvG+nbH1oiI9O1oiK1pbKiLrUk7TOgo5Cdia8ZHR2NrJsbG4vc15fNZTuAwS8rlfTPvtyth3cumtuNUE1ti3q83k32JuM2mcdlWuRKX51Qr+P0IYCGZ84PmzZs3y1VXXTXXuwGAmsP4CAC2qPExLIk3sk40eEaQlXXnnmfddmQ91gwR0xtMGVFBodFmhyIZ27NCqsQIfDJCukIjQCs0HptKqyZzH9YtXRNFK1RKnxzzzPQy47Uy3o+8b51Im95WNN7GjPGaFIzXMzCeQ8pMgzN6YZzktV4n82Sg9f64hltZAV/mZ8r4rMyizgwbcw79cu2zarIDzRzDwUqbZpEDNvdLTm3atEkGBwcn/+zYsWOudwkANYHxEQBsjI8AFpI5v9KczWYlm83O9W4AoOYwPgKAjfERwEIy51eaAQAAAACoVYmvNI+MjMjzzz8/+fdt27bJ1q1bZdGiRbJmzZqKdg4AagnjIwDYGB8B1LLEB82PPPKIvO1tb5v8+2WXXSYiIhdddJHceOONFesYANQaxkcAsFVqfAxCX2RKYJQZ8OXriZRGLpQZ+hWYoV/G9qy+Oa5a4FnhZUZYlJlTZvXPKLQeGxqhX2b/jPCtfGAk5BesoDK941RGrxxgBakVjNegaIWNlW4/rbflp9xC2czth75q8owPkPU+Won7vhWgZb3fxmOtECw7uEs12Y+12lz7Mi9BYI7hYI5tpftwXe3Akvig+ayzznJa7gEAXm8YHwHAxvgIoJZxTzMAAAAAABHmPD0bAJBMIGItwzjJmtama6y5iKVcpxK6FDnszWE7obHGp6qx1rwsERSNqYQlivnx2Jrc6EBsjYjIaCb+n9O0OT90OmO2qBIGeupgqcCPr3HZjstrLWKv/6pqnN7/+CKX7bhMk3XaToX6U8ltlXuNZjP1EAAQrWoHza73fYhELD5frl7K/xA4/oQT5a1nvl3e8MaNcuIbN8ryFStFRGR1d4tZ/28/vF1OPf2Myf/+i4cfFxGRt2w83qk/YRj/w2yy1o//oTeVnx9LVJ8f7UtUP7wv2XIPdXXJPlIuPxCnKuRzyepzyV7PfCGfqL7o8MN8Kt+4TyhK0h8/oXEPTvntJyo371spv/1kb67rj3PrfiFUznEb3iBvOfMdcsIbT5INJ742Ph6+pNms/+6tt8ubT3+rHLn01fHx3od/KyIib3vTsfPTYQCYJ+uP3yCnvfUsOf4NJ8rxJ75Rli1fISIix61aYtZ/899+IG869TQ5Yc2r4+h//teDIiLynjNOm58OA6iY1+WV5kv/z9/Iu959frW7AQALzif+z9/IO89jfASwMIVhOO2KfBgYAU0pK7DI9USuVffqDj926f+Rs9/1bmPb5beowoiMk8Sep++YtJ6HFxh11nQg4xe+2U0rByuwZvy4BWF5RgBX2nhuBStAypj9EhbjT5SnrNCqjJV6puuKxsWM0LzAYb1nVm+swCtjH8b2zEAux7Asa79mlXP41szbzOfrHPDlVmeGiBmf0bDku1H69yRelwfNjz78K3nqySfkN79+VLb++lF58LGnpKGhodrdAoCq+/UjD8nTTz4h//PrR+V/tj4mP3/kSckyPgKA/ObRh+XZp56Ux7f+Wp78zW/kpw8+wvgIvE68Lg+a/7+vfLnaXQCABel6xkcAMP3LV78y+f9TxtVTAAcvvvEAAAAAAETgoBkAAAAAgAivy+nZSb3/D86b9nfX1GwAONhd+Ievjo+p3yWikJoNYC6FMj0YyAr2sVZvsFZ8swKGUmZ8kn7w1FVgDoQSlS4HeNEf/eG0gK93nnryq3VGgpRnBU1ZIWfG8/DS1jUwI6TLesJFvV/fCvgyHmo9j5SVv5XSwWLWShiBEfrlOSwM4hm5Xb4RqpUyQsp8sx/GTs0Pi/GeGf2zPqNmEJgV5mWuGOK4PeN7YIaIzSYIzOqdtQvrdTYf6xgYZjzWCtgr3a/9errhSjMAAAAAABE4aAYAAAAAIALTswFggQnDiPU0k24khjHzL7o2ZnPWtLJSLmdpjeUddY01F6+EtYZjKc9ci3Q6v+h2brmQj69J2Yt6Jua2lVl/ghJtx6nKochlO/YapXO0HYcNuc72q1y/Z/bfAAAzx5VmAAAAAAAicKUZAAAANSGUQMIp8y2ssCMvJrjrACscLDAea4WDTd3vgdkGpUFgkY/VuzVnkJhtRqiUFSImob4u5vlW0JSxE9FBWNZras2esWYwhb7DFCIRkYLxyphhaCVtxuY9YwpIKmMERRmZX0FgNuo263mZSW1Wm/W5cGyzpreYTVY4mK4zg8VmEwRmfRzN/rnuVz/WfgmsUMD4bbmq2kFzGBQldJyqFnpWFF80K+HPhe8yL1DcplhNq09wQT90mC44le8nqy/kJxLVj5uJjNFS1r9AZXjGP3blWP/olRNag1wZxr8n5bdvjbZl692fr8t012n1yZ5q4u0nDRxM+j1x3u6cbBUAAACwMT0bAAAAAIAIr8vp2e84513y6b/+m8m/19fXi4jIj++4b7Lt2n+8Wu6+64757hoAVNXbznmXfPL/vDY+1v1ufPzB7fdOtn3ly1fLvXfdOe99A4BqOvMd58hffvozk38/MD5+7ye3T7ZtufbLcv/dP5v3vgGYW6/Lg+auxYvlpI0nq/apbV1di+ezSwCwIHR1LZYTjfFxatuiru757BIALAiLuhbLG07aqNqnti3q6prPLgGYJ144VzceRhgaGpL29nbpbq9zX34j6T3NCfvkJZyknkq4h1TKfQephJ1JsGkREUknfECGe5rLb39B3dOc9F77ub6nOVm9693KQRDKzt4xGRwclLa2tqQ7WdAOjI8rOhrKjo8uQ6fLJ815CHYodNmWy5hfue1Ual9uY2Cl+u20L6cqlpyqyHZqbMmpIAhl2/6DfXxsmvZdShu/O6zvv9lmfJvsOiNUy3UfZpu1PWOsseqMn8Qp47ea+dwcf4N4xofMfB7Wa5/RHfTSjklORSvkzGEst2qs1yljpZQZYWZF47ejFQ5mhcmY47xbEJiVaGYGaFk/Vp2DuxzHdOP9sYO73FK/7GUgjedrVZl9sXarWwsljw2CUF7qHZnRGFmzV5rvvPsX0r10aVX2PdufPft69srvvfPMivQFAEr95O5fSPeS6oyP5biOnfv27ZX3nfvWOe0LgNenH/30Plm8ZEnkf7fHqcqc8Eq+OV24f1+PXPD776lodwDEq9mD5u6lS2XFipXV7gYAVFwYlj8XHHeeuHvJUllew+NjKPaZ5cn/7nKl3WEZBZffrYHjcgwuk2x8lysmFZpFUKtme6U1yZZctlOpK9bztb95njxYkxYvWVLT4yOA6qjZg+Z9e/dWbd+VuNIMAHOlZ++eandhVvYzRgKYI/t7esr+91q40gxg/tXsQfO573hL5H+r5XuaAWC23vOOt7jd0+u4vfm8X/lgvooKoPre+86zJv//wXJPM4C5V7MHzQAAAHidKZmCbk5It/KZjHsonEORjNs0Uq4BcNaBuVFnt1mBT0ZwlW+FMTkGnxn7tV5V+8SocVBvBaRaAV9Wn633zQzR8uJrrNfJun3BOAlh3uZgJsU6tln3FJlvuFvglR36NfO60HrhzTrHsDHXADIzz2zmoV9mWHBp0yzuYOESKAAAAAAAEThoBgAAAAAgQtWmZ4eSMqes2MVJr6UnvCsu2VK+Tsmt0zafYIHb0EvWGXvds3J9Sbb9op/w/u3E65AmfG+Tbj/p2sKJ520k7Y/765/0Y594XecFtv1qbxcAAACwcKUZAAAAAIAIBIEBAACgJpRO3nNcSj0iJMgxpMsKHTKznYztWcFGRqc9I1nMWv89dJy9aL8u1vM1+mLM6gvNdG89e86z+mcFfFkT74yHmmFoKgjMugboFjwVkYxltFnBZa5Pwngjzewt670wdmHt1XjdzafmPF3PLeDLnD/pvF+3wDA7bMxtH6UzcpPO0J2Kg2YAWGBCibmzYAFOUTd/kJYwE05VTTzzR1QV9yXidtuA05JbLrur0Lpc8728l9Mr6VDk+I7EV7jsq2L9cd3W7HbI7SsAMDeYng0AAAAAQAQOmgEAAAAAiMBBMwAAAAAAEbinGQAAADXBD0tu67ZCpqwb9otWnVvwknn/v1lnBXwZdcbN557xRPyka5xO255jnWugmeMGrbLQN+qs9801REtKg8AcH2f0QyXLif25MAPYfCtEIOE6ttN2YrS5vpFuOVtmEJadpeBYZwSkmcFddjqYw17d66xMh2JJ9xKsAqzM6ErzV7/6VTnkkEOkoaFBTjnlFHnooYdm3gMAOIgwPgKAjfERQK1KfND8ve99Ty677DK54oor5LHHHpMTTjhBzj33XOnp6ZmL/gFAzWB8BAAb4yOAWpb4oPnLX/6y/Pmf/7lcfPHFsn79evna174mTU1N8s1vfnMu+gcANYPxEQBsjI8Aalmie5rz+bw8+uijsmnTpsm2VColZ599tvzyl780H5PL5SSXy03+fXBwUERmt7h0nKR3gMz1WpVO624eqE3YGeu+i/IPSLr9ZPXJX8y57X/SzTutkTkbCT73Sb8iSfu+0Lbv6sDYYd0zU00LbXx0W4O4gtuq1DrFLvty2c487stVpfrEOs2V2ZDT16zG1ml+PY6P1thjfrYd61zb7O1Z9ypbbdZjrZuLq3NPs/Voh1uLI+uM214rfE/zLF478752x3uazRtkX2/3NLveq1yde5pLn+9sxshEB8379+8X3/dl6dKl09qXLl0qTz/9tPmYzZs3y1VXXaXa+4bySXYNANP09vZKe3t7tbsxqZLjYw/jI4BZOJjHx96RiTnpI4DXj5mMkXOenr1p0ya57LLLJv8+MDAga9eule3bty+oAb2coaEhWb16tezYsUPa2tqq3R1ntdjvWuyzCP2eT4ODg7JmzRpZtGhRtbsya4yP1UO/508t9lmkNvvN+Ljw1OLnqBb7LFKb/a7FPovUbr9nM0YmOmhevHixpNNp2bt377T2vXv3yrJly8zHZLNZyWazqr29vb2mXmQRkba2tprrs0ht9rsW+yxCv+dTKrWwlplnfKy9z5AI/Z5PtdhnkdrsN+PjwlOLn6Na7LNIbfa7FvssUrv9nskYmegR9fX1ctJJJ8ndd9892RYEgdx9991y6qmnJt45ABwsGB8BwMb4CKDWJZ6efdlll8lFF10kGzdulJNPPlmuvfZaGR0dlYsvvngu+gcANYPxEQBsjI8Aalnig+YLLrhA9u3bJ5dffrns2bNH3vCGN8gdd9yhwh2iZLNZueKKK8wpNwtVLfZZpDb7XYt9FqHf82kh95nxsXbQ7/lTi30Wqc1+L+Q+vx7HR5Ha7Hct9lmkNvtdi30WeX322wsX2roEAAAAAAAsEAsrKQIAAAAAgAWEg2YAAAAAACJw0AwAAAAAQAQOmgEAAAAAiMBBMwAAAAAAEeb1oPmrX/2qHHLIIdLQ0CCnnHKKPPTQQ/O5+8QeeOABOf/882XFihXieZ788Ic/rHaXYm3evFne9KY3SWtrqyxZskT+4A/+QJ555plqdyvWli1bZMOGDdLW1iZtbW1y6qmnyu23317tbiVy9dVXi+d5cumll1a7K2VdeeWV4nnetD9HH310tbvlZOfOnfKhD31Iurq6pLGxUY4//nh55JFHqt2tiqi18VGEMXK+MD7OH8bHhavWxkjGx/nDGDl/anWMrMT4OG8Hzd/73vfksssukyuuuEIee+wxOeGEE+Tcc8+Vnp6e+epCYqOjo3LCCSfIV7/61Wp3xdn9998vl1xyiTz44INy1113SaFQkHe+850yOjpa7a6VtWrVKrn66qvl0UcflUceeUTe/va3y3vf+1757W9/W+2uOXn44Yfl+uuvlw0bNlS7K06OPfZY2b179+SfX/ziF9XuUqz+/n45/fTTpa6uTm6//XZ58skn5R//8R+ls7Oz2l2btVocH0UYI+cL4+P8YnxceGpxjGR8nD+MkfOr1sbIio2P4Tw5+eSTw0suuWTy777vhytWrAg3b948X12YFREJb7311mp3I7Genp5QRML777+/2l1JrLOzM/yXf/mXancj1vDwcHjkkUeGd911V3jmmWeGn/rUp6rdpbKuuOKK8IQTTqh2NxL77Gc/G77lLW+pdjfmRK2Pj2HIGDnfGB/nBuPjwlTrYyTj4/xjjJwbtThGVmp8nJcrzfl8Xh599FE5++yzJ9tSqZScffbZ8stf/nI+uvC6NTg4KCIiixYtqnJP3Pm+LzfffLOMjo7KqaeeWu3uxLrkkkvkPe95z7TP90L33HPPyYoVK+Swww6TCy+8ULZv317tLsW67bbbZOPGjfLHf/zHsmTJEjnxxBPl61//erW7NWuMj9VVa2Mk4+PcY3xcWBgjq6fWxkcRxsj5UGtjZKXGx3k5aN6/f7/4vi9Lly6d1r506VLZs2fPfHThdSkIArn00kvl9NNPl+OOO67a3Yn1+OOPS0tLi2SzWfnYxz4mt956q6xfv77a3Srr5ptvlscee0w2b95c7a44O+WUU+TGG2+UO+64Q7Zs2SLbtm2TM844Q4aHh6vdtbJefPFF2bJlixx55JFy5513ysc//nH55Cc/Kd/+9rer3bVZYXysnloaIxkf5wfj48LDGFkdtTQ+ijBGzpdaHCMrNT5m5qh/WAAuueQSeeKJJxb8vQYHrFu3TrZu3SqDg4Nyyy23yEUXXST333//gh30duzYIZ/61KfkrrvukoaGhmp3x9l55503+f83bNggp5xyiqxdu1b+7d/+TT760Y9WsWflBUEgGzdulC984QsiInLiiSfKE088IV/72tfkoosuqnLvUItqaYxkfJwfjI/Aq2ppfBRhjJwvtThGVmp8nJcrzYsXL5Z0Oi179+6d1r53715ZtmzZfHThdecTn/iE/OQnP5F7771XVq1aVe3uOKmvr5cjjjhCTjrpJNm8ebOccMIJ8k//9E/V7lakRx99VHp6euSNb3yjZDIZyWQycv/998t1110nmUxGfN+vdheddHR0yFFHHSXPP/98tbtS1vLly9U/fsccc8yCnxYUh/GxOmptjGR8rA7Gx+pjjJx/tTY+ijBGVkstjJGVGh/n5aC5vr5eTjrpJLn77rsn24IgkLvvvrsm7jeoJWEYyic+8Qm59dZb5Z577pFDDz202l2asSAIJJfLVbsbkd7xjnfI448/Llu3bp38s3HjRrnwwgtl69atkk6nq91FJyMjI/LCCy/I8uXLq92Vsk4//XS19MWzzz4ra9eurVKPKoPxcX4dLGMk4+P8YHysPsbI+XOwjI8ijJHzpRbGyEqNj/M2Pfuyyy6Tiy66SDZu3Cgnn3yyXHvttTI6OioXX3zxfHUhsZGRkWlnTrZt2yZbt26VRYsWyZo1a6rYs2iXXHKJ3HTTTfKjH/1IWltbJ+/3aW9vl8bGxir3LtqmTZvkvPPOkzVr1sjw8LDcdNNNct9998mdd95Z7a5Fam1tVff5NDc3S1dX14K+/+czn/mMnH/++bJ27VrZtWuXXHHFFZJOp+UDH/hAtbtW1qc//Wk57bTT5Atf+IK8//3vl4ceekhuuOEGueGGG6rdtVmrxfFRhDFyvjA+zh/Gx4WpFsdIxsf5wxg5f2pxjKzY+Djr/O0EvvKVr4Rr1qwJ6+vrw5NPPjl88MEH53P3id17772hiKg/F110UbW7Fsnqr4iE3/rWt6rdtbI+8pGPhGvXrg3r6+vD7u7u8B3veEf405/+tNrdSqwWlgu44IILwuXLl4f19fXhypUrwwsuuCB8/vnnq90tJz/+8Y/D4447Lsxms+HRRx8d3nDDDdXuUsXU2vgYhoyR84Xxcf4wPi5ctTZGMj7OH8bI+VOrY2QlxkcvDMNwRoftAAAAAAAc5OblnmYAAAAAAGoRB80AAAAAAETgoBkAAAAAgAgcNAMAAAAAEIGDZgAAAAAAInDQDAAAAABABA6aAQAAAACIwEEzAAAAAAAROGgGAAAAACACB80AAAAAAETgoBkAAAAAgAj/PxGm7ufUmvFHAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1200x400 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "numpix = 64  # number of pixels per axis of the image to be modelled\n",
    "\n",
    "# here we define the numerical options used in the ImSim module. \n",
    "# Have a look at the ImageNumerics class for detailed descriptions.\n",
    "# If not further specified, the default settings are used.\n",
    "kwargs_numerics = {'point_source_supersampling_factor': 1}\n",
    "\n",
    "\n",
    "def simulate_rgb(ConfigList, size, kwargs_numerics):\n",
    "    band_b, band_g, band_r = ConfigList\n",
    "    kwargs_b_band = band_b.kwargs_single_band()\n",
    "    kwargs_g_band = band_g.kwargs_single_band()\n",
    "    kwargs_r_band = band_r.kwargs_single_band()\n",
    "    \n",
    "    # set number of pixels from pixel scale\n",
    "    pixel_scale = kwargs_g_band['pixel_scale']\n",
    "    numpix = int(round(size / pixel_scale))\n",
    "\n",
    "    sim_b = SimAPI(numpix=numpix, kwargs_single_band=kwargs_b_band, kwargs_model=kwargs_model)\n",
    "    sim_g = SimAPI(numpix=numpix, kwargs_single_band=kwargs_g_band, kwargs_model=kwargs_model)\n",
    "    sim_r = SimAPI(numpix=numpix, kwargs_single_band=kwargs_r_band, kwargs_model=kwargs_model)\n",
    "\n",
    "    # return the ImSim instance. With this class instance, you can compute all the\n",
    "    # modelling accessible of the core modules. See class documentation and other notebooks.\n",
    "    imSim_b = sim_b.image_model_class(kwargs_numerics)\n",
    "    imSim_g = sim_g.image_model_class(kwargs_numerics)\n",
    "    imSim_r = sim_r.image_model_class(kwargs_numerics)\n",
    "\n",
    "\n",
    "    # turn magnitude kwargs into lenstronomy kwargs\n",
    "    kwargs_lens_light_g, kwargs_source_g, kwargs_ps_g = sim_b.magnitude2amplitude(kwargs_lens_light_mag_g, kwargs_source_mag_g, kwargs_ps_mag_g)\n",
    "    kwargs_lens_light_r, kwargs_source_r, kwargs_ps_r = sim_g.magnitude2amplitude(kwargs_lens_light_mag_r, kwargs_source_mag_r, kwargs_ps_mag_r)\n",
    "    kwargs_lens_light_i, kwargs_source_i, kwargs_ps_i = sim_r.magnitude2amplitude(kwargs_lens_light_mag_i, kwargs_source_mag_i, kwargs_ps_mag_i)\n",
    "\n",
    "\n",
    "    image_b = imSim_b.image(kwargs_lens, kwargs_source_g, kwargs_lens_light_g, kwargs_ps_g)\n",
    "    image_g = imSim_g.image(kwargs_lens, kwargs_source_r, kwargs_lens_light_r, kwargs_ps_r)\n",
    "    image_r = imSim_r.image(kwargs_lens, kwargs_source_i, kwargs_lens_light_i, kwargs_ps_i)\n",
    "\n",
    "    # add noise\n",
    "    image_b += sim_b.noise_for_model(model=image_b)\n",
    "    image_g += sim_g.noise_for_model(model=image_g)\n",
    "    image_r += sim_r.noise_for_model(model=image_r)\n",
    "\n",
    "    # and plot it\n",
    "\n",
    "    img = np.zeros((image_g.shape[0], image_g.shape[1], 3), dtype=float)\n",
    "    #scale_max=10000\n",
    "    def _scale_max(image): \n",
    "        flat=image.flatten()\n",
    "        flat.sort()\n",
    "        scale_max = flat[int(len(flat)*0.95)]\n",
    "        return scale_max\n",
    "    img[:,:,0] = plot_util.sqrt(image_b, scale_min=0, scale_max=_scale_max(image_b))\n",
    "    img[:,:,1] = plot_util.sqrt(image_g, scale_min=0, scale_max=_scale_max(image_g))\n",
    "    img[:,:,2] = plot_util.sqrt(image_r, scale_min=0, scale_max=_scale_max(image_r))\n",
    "    data_class = sim_b.data_class\n",
    "    return img, data_class\n",
    "\n",
    "size = 6. # width of the image in units of arc seconds\n",
    "\n",
    "img_des, coords_des = simulate_rgb(des, size=size, kwargs_numerics=kwargs_numerics)\n",
    "img_lsst, coords_lss = simulate_rgb(lsst, size=size, kwargs_numerics=kwargs_numerics)\n",
    "img_roman, coords_roman = simulate_rgb(roman, size=size, kwargs_numerics=kwargs_numerics)\n",
    "\n",
    "from lenstronomy.Plots.plot_util import coordinate_arrows, scale_bar\n",
    "\n",
    "f, axes = plt.subplots(1, 3, figsize=(12, 4))\n",
    "\n",
    "ax = axes[0]\n",
    "ax.imshow(img_des, aspect='equal', origin='lower', extent=[0, size, 0, size])\n",
    "ax.set_title('DES')\n",
    "# coordinate_arrows(ax, d=1, coords=coords_des, color='w', font_size=15, arrow_size=0.05)\n",
    "scale_bar(ax, d=size, dist=1., text='1\"', color='w', font_size=15, flipped=False)\n",
    "\n",
    "ax = axes[1]\n",
    "ax.imshow(img_lsst, aspect='equal', origin='lower', extent=[0, size, 0, size])\n",
    "ax.set_title('LSST')\n",
    "scale_bar(ax, d=size, dist=1., text='1\"', color='w', font_size=15, flipped=False)\n",
    "\n",
    "ax = axes[2]\n",
    "ax.imshow(img_roman, aspect='equal', origin='lower', extent=[0, size, 0, size])\n",
    "ax.set_title('Roman')\n",
    "scale_bar(ax, d=size, dist=1., text='1\"', color='w', font_size=15, flipped=False)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "tf213",
   "language": "python",
   "name": "tf213"
  },
  "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}