{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Point plus extended lens simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook shows how to use SLSim to simulate point plus extended source lens. For \n",
    "\n",
    "this, it uses a simple point plus extended source catalog. Catalog has skypy galaxaly \n",
    "\n",
    "catalog as host galaxies and fake quasars as point source."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulate source and deflector population"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from astropy.cosmology import FlatLambdaCDM\n",
    "from astropy.units import Quantity\n",
    "from slsim.lens_pop import LensPop\n",
    "from slsim.image_simulation import lens_image\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define a cosmology\n",
    "cosmo = FlatLambdaCDM(H0=70, Om0=0.3)\n",
    "\n",
    "\n",
    "# define a sky area\n",
    "sky_area = Quantity(value=1, unit=\"deg2\")\n",
    "\n",
    "\n",
    "# define limits in the intrinsic deflector and source population (in addition to the\n",
    "# skypy config\n",
    "# file)\n",
    "kwargs_deflector_cut = {\"band\": \"g\", \"band_max\": 23, \"z_min\": 0.01, \"z_max\": 2.5}\n",
    "kwargs_source_cut = {\"band\": \"g\", \"band_max\": 26, \"z_min\": 0.1, \"z_max\": 5.0}\n",
    "\n",
    "## create a point plus extended source lens population.\n",
    "pes_lens_pop = LensPop(\n",
    "    deflector_type=\"all-galaxies\",\n",
    "    source_type=\"quasar_plus_galaxies\",\n",
    "    kwargs_deflector_cut=kwargs_deflector_cut,\n",
    "    kwargs_source_cut=kwargs_source_cut,\n",
    "    variability_model=\"sinusoidal\",\n",
    "    kwargs_variability={\"amp\", \"freq\"},\n",
    "    kwargs_mass2light=None,\n",
    "    skypy_config=None,\n",
    "    sky_area=sky_area,\n",
    "    cosmo=cosmo,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Select a lens at random"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define a lens cut\n",
    "kwargs_lens_cut = {\"min_image_separation\": 0.8, \"max_image_separation\": 10}\n",
    "# Select a lens at random\n",
    "pes_lens_class = pes_lens_pop.select_lens_at_random(**kwargs_lens_cut)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define a psf kernel and tranform matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load your psf kernel and transform matrix\n",
    "psf_kernel = np.load(\"../..//tests/TestData/psf_kernels_for_deflector.npy\")\n",
    "transform_matrix = np.array([[0.2, 0], [0, 0.2]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulate lens image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/sanchez/.virtualenvs/slsim/lib/python3.12/site-packages/lenstronomy/Data/psf.py:80: UserWarning: Input PSF model has at least one negative element, which is unphysical except for a PSF of an interferometric array.\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "source": [
    "# Simulate a lens image\n",
    "image_lens = lens_image(\n",
    "    lens_class=pes_lens_class,\n",
    "    band=\"i\",\n",
    "    mag_zero_point=27,\n",
    "    num_pix=64,\n",
    "    psf_kernel=psf_kernel,\n",
    "    transform_pix2angle=transform_matrix,\n",
    "    exposure_time=30,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x14fcc4320>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAGfCAYAAAAZGgYhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAzU0lEQVR4nO3df3Bc1X338c+uJa1k2Vr5pyTHsus8MbGB2IABo5q0FJR4/GQYU9yUZMjUTZkwUOGATSdBnQBpJ0GETIJDYkyg1CRPQ524M4aQFlzGxOIhsQ0WMAGcOIa4tYItmV9aybK1krX3+YOyD2LPcXysszq7q/drZmfg7vXZc+69q6+u79ffbyyKokgAAIyxeOgJAADGJwIQACAIAhAAIAgCEAAgCAIQACAIAhAAIAgCEAAgCAIQACAIAhAAIAgCEAAgiLJ8DbxhwwZ985vfVFdXlxYvXqzvfve7uvDCC//gn8tkMjp06JAmT56sWCyWr+kBAPIkiiL19fVp1qxZisdPcp8T5cHmzZujioqK6J//+Z+jV155JfrCF74Q1dbWRt3d3X/wz3Z2dkaSePHixYtXkb86OztP+vM+FkX+i5EuXbpUF1xwgb73ve9JeveuprGxUWvWrNEtt9xy0j+bSqVUW1uri/W/VaZy31PLL9sdWyHVezXNsZDm50s+1zkejmExXMuFZDxcEw5OaEjP6D/U09OjZDJp3c/7X8ENDg6qo6NDra2t2W3xeFzNzc3auXNnzv7pdFrpdDr7/319ff8zsXKVxUokAKmALkTjHAtofr7kc53j4RgWw7VcSMbDNeHif5b+hx6jeE9CePPNNzU8PKy6uroR2+vq6tTV1ZWzf1tbm5LJZPbV2Njoe0oAgAIUPAuutbVVqVQq++rs7Aw9JQDAGPD+V3DTp0/XhAkT1N3dPWJ7d3e36uvrc/ZPJBJKJBK+pxFGMfydr485+ng+kO9nDPk8FyHOs0tGqI/5FcO17KKYr7cS5v0OqKKiQkuWLNH27duz2zKZjLZv366mpibfHwcAKFJ5+XdA69at0+rVq3X++efrwgsv1Pr169Xf36/Pf/7z+fg4AEARyksAuuqqq/TGG2/otttuU1dXl8455xw98cQTOYkJAIDxKy//Dmg0ent7lUwmdYlWFl8a9nhRDM+ASs1YPwMqNVxvY+pENKQdelSpVEo1NTXW/YJnwQEAxqe81YJDCchnLb5i+M2z0H9rLpR5FINiPlYl/DcO3AEBAIIgAAEAgiAAAQCCIAABAIIgCQF2xfzg1odCesibx3MRK6/I/bgTQ2M+D1g4HPNYmflHenTihK/ZeMUdEAAgCAIQACAIAhAAIAgCEAAgCAIQACAIsuCQf6YMsXxnU431Z/oaO0DZlWho8NTH9vSZ40KAY2LNdnOdyxh9f7gDAgAEQQACAARBAAIABEEAAgAEQQACAARBFtx4EyIjLUQmlOkzCylTq5CajIVo9x3iOhxrhbQe17mM0dy5AwIABEEAAgAEQQACAARBAAIABEEAAgAEQRZcsXPNhCqUzBzHeXvp9Fgoa5cUTySM2zMDA6c+SBF0W3U6z4WUpYgxwR0QACAIAhAAIAgCEAAgCAIQACAIkhAKkUuZkmJ9QOs4b6dkAxsfTblOtr8DL+spJD7K//gqIVSs34liXY9x3jHpFKbNHRAAIAgCEAAgCAIQACAIAhAAIAgCEAAgCLLgClFRZr3ILZss35lnLk3WbPK4nrxm9dkUSskd17kUa3aYq2Jdj2nep7gW7oAAAEEQgAAAQRCAAABBEIAAAEEQgAAAQZAFB3euGWz5HDtEjTyXZmquXNZZSFlTHs5brLzCvOvQ4KjHtgpQH9DaXHF4eNRju/LS6HEUuAMCAARBAAIABEEAAgAEQQACAARBAAIABOEcgJ5++mldfvnlmjVrlmKxmB555JER70dRpNtuu00NDQ2qqqpSc3Oz9u/f72u+41csZn752t/HGFF06i9XvsbxwXRMXNZ+snn7WKePcx9g7OjEkPHlhW3etuPtur/TOk8YXyGucetcxohzAOrv79fixYu1YcMG4/t33XWX7rnnHt13333avXu3qqurtXz5cg0MDIx6sgCA0uH874BWrFihFStWGN+Lokjr16/XV77yFa1cuVKS9MMf/lB1dXV65JFH9JnPfCbnz6TTaaXT6ez/9/b2uk4JAFCEvD4DOnDggLq6utTc3JzdlkwmtXTpUu3cudP4Z9ra2pRMJrOvxsZGn1MCABQorwGoq6tLklRXVzdie11dXfa9D2ptbVUqlcq+Ojs7fU4JAFCggpfiSSQSSiQSoacBABhjXgNQfX29JKm7u1sNDQ3Z7d3d3TrnnHN8ftTYcel+WUifWegdJ31183QZx9d6THXMfNXUCrEe0zR8rcdHV1lXIb6zLnxlJLrUJCyk9b+P17+Cmzdvnurr67V9+/bstt7eXu3evVtNTU0+PwoAUOSc74COHj2qV199Nfv/Bw4c0IsvvqipU6dqzpw5uummm/S1r31N8+fP17x583Trrbdq1qxZuuKKK3zOGwBQ5JwD0J49e/Rnf/Zn2f9ft26dJGn16tV66KGH9KUvfUn9/f269tpr1dPTo4svvlhPPPGEKisr/c0aAFD0YlFUWH852Nvbq2QyqUu0UmWx8tDTKd5nQKP9vFCfaVNIz4BM0yikZ0D57Fnj61/J+7jG83nu8/mdGAfPgE5EQ9qhR5VKpVRTU2PdL3gWXMErpAvOh0J5QFtIPxBsH+nQOMzaTMwXH03wHI5h3tfjQz4TbXwE93wGGp/jB0QxUgBAEAQgAEAQBCAAQBAEIABAEAQgAEAQxZUFVyAphs4KaX55zBzKe+quCw8putZ5e8gkdD5W+UxbNmyPTZhg3NXbufSQNeY0R0/fQdtnnvI8JH8/xwrln0mMYh7cAQEAgiAAAQCCIAABAIIgAAEAgiAAAQCCKK4suGItyFlIhQ1dPtOxBlVes91cM+8Mtcy8ZXZ5KALpXGstnzXvDMcl75mLHuqY2eZouyZcxvC1v3mQAsqKtclnbbv34Q4IABAEAQgAEAQBCAAQBAEIABAEAQgAEERxZcHlUz4zU0KM7ZKx4qtbpA+u2WQO9bZcM89cMuycecjq89Up05gxGKKun6/1GOYYr6w85X1PymEu3rIu89n51OE6zMe55w4IABAEAQgAEAQBCAAQBAEIABBE4SYhxGK5D9+KoYRFofBRciefXBMf8llyyDaELdnAMLavh/a2/U0P0TPptNPYTvOwrD1WXmHef2jQaXwv5X8cronMwIDb2DYO11Vem/flmZdEm1PAHRAAIAgCEAAgCAIQACAIAhAAIAgCEAAgiMLNgosiSR/I/shnYzcfCml+IeYSoGSIlWkurmO4ZDzZsoY8nQcvWVwuc7E2HRwa/TxUQOV/8inf30Ef17iNh6aLp4I7IABAEAQgAEAQBCAAQBAEIABAEAQgAEAQhZsFZ1IojdACZJg5Zwj5mIvjOq0NuDw0U3Ma21E+s69iZeWWsS3ZZB5q4dmOlU0+6365NNhznofD9encXNDH9yef2W628QspE/cUcAcEAAiCAAQACIIABAAIggAEAAiCAAQACKK4suDyqcA7iOazW6Trepyzxkxzcawb55ytZBjfdd62/eOTqnO2Dad6LWObs92s2XGWzqLxiRNztmWOHTOPYTlW8UTCuF2xU/89NDbBvK/z9ekjU83l41yuzZNtD5FNls+fTS7rzMPauQMCAARBAAIABEEAAgAEQQACAAThFIDa2tp0wQUXaPLkyZo5c6auuOIK7du3b8Q+AwMDamlp0bRp0zRp0iStWrVK3d3dXicNACh+Tllw7e3tamlp0QUXXKATJ07o7//+7/XJT35Se/fuVXX1u5lBa9eu1b//+79ry5YtSiaTuuGGG3TllVfqF7/4RV4WEMQ46XRorb/mmlHkMnfHubhkcNkz0tySQYd7UoZBLPO2ZLtZ2TqRptOGjZaafOUV5jGGM+b9DZltGdPnSYpsp95yflyOret15aU+oIfae1auYxR6hl0exKLo9GfwxhtvaObMmWpvb9ef/MmfKJVKacaMGXr44Yf1F3/xF5Kk3/zmN1q4cKF27typiy666A+O2dvbq2QyqUu0UmUxxy/vWCmGAOShgOG4CUCuxTtN6/cUgFzmaE0ftwQgG5cAZB/DLQCZ5l4UAcjDtew8lyJ0IhrSDj2qVCqlmpoa636jegaUSr37m+DUqVMlSR0dHRoaGlJzc3N2nwULFmjOnDnauXOncYx0Oq3e3t4RLwBA6TvtAJTJZHTTTTdp2bJlOvvssyVJXV1dqqioUG1t7Yh96+rq1NXVZRynra1NyWQy+2psbDzdKQEAishpB6CWlha9/PLL2rx586gm0NraqlQqlX11dnaOajwAQHE4rVI8N9xwg372s5/p6aef1uzZs7Pb6+vrNTg4qJ6enhF3Qd3d3aqvrzeOlUgklLCVCMkHH3+36+t5jMv4rvP20pDOU5a+j7JAjnOJVeQ+e7GVkbE977A9v4lXVuZssz2/sD3TsZXFiVVOMm83Pac52u/0mTamxALb/OzPaSwleiyJD06leyzXhJdGeq7fKw+la/LZANFZ4JJDTt/qKIp0ww03aOvWrXrqqac0b968Ee8vWbJE5eXl2r59e3bbvn37dPDgQTU1NfmZMQCgJDjdAbW0tOjhhx/Wo48+qsmTJ2ef6ySTSVVVVSmZTOqaa67RunXrNHXqVNXU1GjNmjVqamo6pQw4AMD44RSANm7cKEm65JJLRmzftGmT/vqv/1qSdPfddysej2vVqlVKp9Navny57r33Xi+TBQCUDqcAdCr/ZKiyslIbNmzQhg0bTntSAIDSRy04AEAQNKT7Q3yU5PDBV9M4U+aQLcvI0hzNxtQ0TTKXkXFtAqfInE1lEw0aMsEsY9gyvqzZcRW564yVWyoeHD9u3m7bf8icwRa5JbYZWasyGI6LLXvNynJsrc3xDJmEmYEB476upZJM+9vGduYhO8z12neqEOE6P5dsPxrSAQBKBQEIABAEAQgAEAQBCAAQBAEIABDE+MuC85El4qt+Uh6bXnmpK+W4zky/uTaZ0zptNcUsGXmmbCpJxiyzaMCS1VZpzoJz+u3Mkr1mzeDKWLLMLNlxsYlVOduit3uM+05IVhu3D6dOvdVJPGGedybtVtvOie1adqz55nTt57MWnOPYXmrb+VKIteAAAPCFAAQACIIABAAIggAEAAiCAAQACKK4suDy2b3PQ6dDL2PLnDnlnNXmIVvHpXaYZM/iMY1j69pp61oanTDP0VavLT5hwqmPbat7ZstIM4xjrXlWmzR/puM1lDFkvNmy92zriU9y6Lbab65hFzMc15N9pi2r0XTe7GPnMTvMU7aXl++srzpuBgXVhfV9uAMCAARBAAIABEEAAgAEQQACAARBAAIABFFcWXD5rE9UQGN7yUzx0rnRsQ2nrcaVYRxbhp01g8vWbdXU+dS23dYRNVlj3K7J5ppqevPt3DEs2W6yZHbJ1m3Vkq0UT07O3bcqtz6cJEWDlk62tmwy0xwtWXC2rDbXjrUumZEu2aL2DzTP25YZae2g6lLHLZ9Zu45C1IY8FdwBAQCCIAABAIIgAAEAgiAAAQCCKNwkhFgs96GXy8Mu12ZvLg2oHB+6Wctg2B4KeyiXYy11Y9rfVlrH8uAyVl7h9JnGfW0N5qotD/5tZX4s5WhMDeJik3If5EuyJwTYGqRNn5q78egx89hxy8PvyeayONFES4M9QwO76O2UeV9HxjI/FZbrynKsbHyUgLFdb9bEB0PCge16s5V4cp23qYyQawkh52NluD6t5YwKJLHpg7gDAgAEQQACAARBAAIABEEAAgAEQQACAARRuFlwUSRpFFkXeW4a5zSVPDZ9sjd2s2XDmErUeGrK5dBQzJpJZ8uyspZSMWdCxWzldUz7WrKPhmfPMO8/mHs+BxunmPfNmI9tWY+51MvxD5mz4yrfyN0/Vm4+3vEj7xi3q9pczih23DAXS/ZefJohA1BSlOo1b3e49uMJc0ajremg0zVu41JaR7L+nHD6jvsYQzLOPa8N5ijFAwAoFQQgAEAQBCAAQBAEIABAEAQgAEAQhZsF58JDvTYvdeZ81UrysB6n+lGOtd1cM+9cmo/ZattZ2WrEGdaZ6bfUa6ubbtwc7zdnX2Wqc7O1bNlubywyN41LHjCvc6ja/DvhYE1ujbzq183zq0j1G7crba6HFmuYmbMt6n7TuG90zNyoztZI0Ja9aKrN5tRgLt9cakPattvGCNCQzsrlZxm14AAApYIABAAIggAEAAiCAAQACIIABAAIooDSTkYhn1klrp1VfYzhYT3WzDZD9pFT7Sy5d100zcU12y0+tdb8hqVOWDQltxZcbJK5Flo00VyDbGiaeX+T1z9u7mQ6XGU+l6mPmbfX/sr8O2HVG7lZZkOTzV/feIO5Ll3ZYUuNuB5zHTeT2ATL76yuHVQN3wnr9ePY9fdUP+/dQcznwbnrr4+fQfnOrnUZ20dm8SngDggAEAQBCAAQBAEIABAEAQgAEIRTEsLGjRu1ceNG/dd//Zck6ayzztJtt92mFStWSJIGBgZ08803a/PmzUqn01q+fLnuvfde1dXVeZ/4aXN90OfjwZuv5ngOY1tL3RgerpoSE07K1hzO4cFtfKLlAb+lEZqptI4kqdL88F+9ueVoIksSwvAkSxJCtfnr0d+Qm4QxwfKsfcI55gf8J3rN8+5ZYh7oxCu5c0weMJ/78oQlScSSKBCLG5rgRX3GfTO9R81j2JIQbKWSDIkFtgQHa3KCSzKMYyKQa2KOF4VUomeMON0BzZ49W3feeac6Ojq0Z88eXXrppVq5cqVeeeUVSdLatWv12GOPacuWLWpvb9ehQ4d05ZVX5mXiAIDi5nQHdPnll4/4/69//evauHGjdu3apdmzZ+vBBx/Uww8/rEsvvVSStGnTJi1cuFC7du3SRRdd5G/WAICid9rPgIaHh7V582b19/erqalJHR0dGhoaUnNzc3afBQsWaM6cOdq5c6d1nHQ6rd7e3hEvAEDpcw5AL730kiZNmqREIqHrrrtOW7du1Zlnnqmuri5VVFSotrZ2xP51dXXq6uqyjtfW1qZkMpl9NTY2Oi8CAFB8nAPQRz/6Ub344ovavXu3rr/+eq1evVp79+497Qm0trYqlUplX52dnac9FgCgeDiX4qmoqNBHPvIRSdKSJUv03HPP6Tvf+Y6uuuoqDQ4OqqenZ8RdUHd3t+rr663jJRIJJRLmLKRRCVDWwtZQy9ocziaPzfFM2T3O847MzcdsGU9xw/m1lWiJVZkbuFlZMqEyvblZXNGsacZ9I0v21YQB83qGE7nHq3+BeT1rzjD/9fOWznON2/vT5kzCvv+Vu86yY+bzNrHL8Ro/PpCzKRq0ZFFast2sjQFtTQoNx9x2vblk0r37hkMZGV9lsnyUrglQsivI2O8z6n8HlMlklE6ntWTJEpWXl2v79u3Z9/bt26eDBw+qqalptB8DACgxTndAra2tWrFihebMmaO+vj49/PDD2rFjh7Zt26ZkMqlrrrlG69at09SpU1VTU6M1a9aoqamJDDgAQA6nAHTkyBH91V/9lQ4fPqxkMqlFixZp27Zt+sQnPiFJuvvuuxWPx7Vq1aoR/xAVAIAPcgpADz744Enfr6ys1IYNG7Rhw4ZRTQoAUPqoBQcACKJgG9LFysoUi42cnlM2WZD6a47Zbj6a3TkyZg5ZstpsrFlJDrW5rGNkLBl2tvNpy6arm567rfe4cd9MebVx+3ClOctqwoBhLkPm3+W+/8rFxu11teZaa31dk43bZfjIij7zMclUWOrpDZhr/pnOj2vjueG33jZ/pqU+YMZw3qwZcx4a0lkzPS3XoX0gD9m1hZTtZmE6XtZjNYr5cQcEAAiCAAQACIIABAAIggAEAAiCAAQACKJgs+CiEycUjSZLzFfnU5c55DO7xVdtO0PGmzW7xfKZ1hpcFqb6YTFbFly5JVspZW7TEZvdYP7Qd1K526YkjbvGh8yZd5Vv5tZIk6ThRG69uoo3zPOO/958rLoT5sy7cnMimKq6c8/FxCPmrLaKI7ndYCUpqjF/ZiyeO3Z01DyG7ZqIV5vHtnK5DiNLR1RbB15Dh19rhqrtGnfpturK8WeNlxqTrjUjXbMDTxN3QACAIAhAAIAgCEAAgCAIQACAIAhAAIAgCjYLbtR81U9yGcdX/TkPc7dlCLnUfXPpZim5ZeVEw+Z52DJ+YtWWuRw314KLZhq6n56wdGztMWd8DcydYh7bsPza3xh3VbrWvD1Tbj73VW+Yz335sdy5lxm2vTu45Ry/Ya7XlkkbsulsY1jYsheHewzZiDKfZ9fMM1O227tzyb32bfvauGbNGXe1rcexppqXjDTnbNnRZujGjPULP4g7IABAEAQgAEAQBCAAQBAEIABAEKWbhOCLh0Z1IRrP2fh4oGlqJiZJ8UTC/AcMn2lNZDhuLn8Tnz7VvL+lZEwsnjt+lDCXxRlsNCcbVKQsD7kNp/n4DPPY0/aaxzhRZUmqGDZfQ5Vdx4zbjWMcNTfei2przPsbyhZFA5ZGf7ZSSUNuJWqcrkPXMjKmhAPXhB/XMlSm0kI+yvZIfhKq8pjwZBzjFMflDggAEAQBCAAQBAEIABAEAQgAEAQBCAAQxPjLgvPVqM6HPDbBi04MmXc3ZfFYyvPYsnhsZX5s2XGmzzQ1qZOkzDFztlfsqCXDzrL+zJRJp7zv8enmucSnmLe/Mz83g63MkqQ2sdt2js2bE++Yj+GEnqM526JKS7klQ4M5SYq63zRPxZSRaMlSs59jy48SD00NXcvoOJXLsTV7s6zf9r3K688JHxlsPsqE5WGN3AEBAIIgAAEAgiAAAQCCIAABAIIgAAEAghh/WXCWTA5bZpdTBo6vekujbgZlH8OUxWOr4WbNBLIck/jkycbtmaO5GVy2WnC2Jl6uMhNzz2d6qnmdx6eb59JzlqUp26TcTLBoyDxGprzSuL3hl7nHRJJiacsxNx2vt83N3qzZi5bMQ2PGm+082BoJVpqPbdzWlM2QTWedtzXDzlJP0OE769oczilrLkRTTF/G6DO5AwIABEEAAgAEQQACAARBAAIABEEAAgAEMf6y4Gy101zrTY0xa/aNY+aQKVvHpYbbyZiy3WzjxKonmgcZNNfaivotNeKS5i6fE94yzGWKJQuuznxNrF/+f4zb55S9k7Pt5tc+bdz3QGK6cXv6N+a5VP/6iHG76biYMskkKVZdbdyeeett43bTNRS3HNfIcn4yls601kw1h46otn1jZQ6/P3vqiGrtcmra37UTcohst8C4AwIABEEAAgAEQQACAARBAAIABDH+khDy2Zgpjw8RXR7ans7+xjFszbcsXJIWht/OfZAv2UvxxGuT5v0nVhm3R2W54xyrM1/u6enmY3UsY04UOMdQuuis2sPGfQ89M9u4fbDGXNKmqsacnBHresu43WS425zIEKsyH6tYee55syV9uF4Tihwe2jtyShxy/G66loSyJicYB3dcu8v+RZbIwB0QACAIAhAAIAgCEAAgCAIQACAIAhAAIIhRZcHdeeedam1t1Y033qj169dLkgYGBnTzzTdr8+bNSqfTWr58ue69917V1dX5mO/Y85FV4qMMiK9GeoaxbRk/PjLpJHPmlDXLyFa6JW1eT6za0iDt2EDOturD5gyuviPmY/iztxYbt88o+78527a9ttC4b5Qwn7fKNy0lh8otWYCG85axNYezNRg8fty43XbMjbs6XivxSnNDPlPWmL3kjjm7Mp/ls1wb1Tl9r1wy5jyOM2q+Gm6+z2nfAT333HP6/ve/r0WLFo3YvnbtWj322GPasmWL2tvbdejQIV155ZWnPUEAQGk6rQB09OhRXX311XrggQc0ZcqU7PZUKqUHH3xQ3/72t3XppZdqyZIl2rRpk375y19q165d3iYNACh+pxWAWlpa9KlPfUrNzc0jtnd0dGhoaGjE9gULFmjOnDnauXOncax0Oq3e3t4RLwBA6XN+BrR582Y9//zzeu6553Le6+rqUkVFhWpra0dsr6urU1dXl3G8trY2/cM//IPrNAAARc7pDqizs1M33nijfvSjH6nS8nDRVWtrq1KpVPbV2dnpZVwAQGFzugPq6OjQkSNHdN5552W3DQ8P6+mnn9b3vvc9bdu2TYODg+rp6RlxF9Td3a36+nrjmIlEQglT1k4slpt14SFrLAjHplcuc7fW5nKoH5XPbDcba6aWpcmaMuaMr+HXzTXYJjR+KHfboHmM+t3mbKpfTltg3P6L8jNyts2Ya65tV/68uf5a2TFzBtOE379h3K6q3F/4YgO5mX6SlBmwNRi0fN0d6p5l+s2N52xj25odmq5xe9NF83Xl2qTRSQH9/BjzbDebPBwTpwB02WWX6aWXXhqx7fOf/7wWLFigL3/5y2psbFR5ebm2b9+uVatWSZL27dungwcPqqmpyd+sAQBFzykATZ48WWefffaIbdXV1Zo2bVp2+zXXXKN169Zp6tSpqqmp0Zo1a9TU1KSLLrrI36wBAEXPezuGu+++W/F4XKtWrRrxD1EBAHi/UQegHTt2jPj/yspKbdiwQRs2bBjt0ACAEkYtOABAEIXbETWKJJ1i1oVDdoY1c8a19pNxcMestnxm2uSxZpVr9pFpf1s2le0Yxiznx3reBnMzpxK/M2eYDX1oqnH7R35szr4aTObWjhuuMo8x6dWUcXvsmDk7LLKct+jNt3O32Wq7WUS2ZEeXumyW2oM28Ymnnu3omu1l3d9HVqyH77K37DU6ogIA4BcBCAAQBAEIABAEAQgAEAQBCAAQROFmwZl4qZ2Wx7pKvjJQxri2nWstOOcMqchcg804hrXjpnmM+MSJ5v0NWXax6mrjvmVHzC1Aompzwd3Kody5TDjSY9xXGct5K7Os01I7LVaR2xU009dnHtv6PbFkKRrOWzxhyxY1nwdbvTZrfUAf2aU2Pr4rrvUb8zXGycbxIQ9dTl1wBwQACIIABAAIggAEAAiCAAQACKK4khAKpMyEt0ZYLg8jbSVqHMvo+DiGts90GsPWHM3C+pB7yNxMzvRgPTZkeQiftnT37X7TuDk+Iff3tsjWYC9lSXCwlpGx/E5oSOSwHkPLGDHDvG1zyaRPvbngu4N4eODueo3bkmdCNKh0Gb9Afo5JCj4X7oAAAEEQgAAAQRCAAABBEIAAAEEQgAAAQRRXFpyNS9ZLIZXzsXymKbvJNWvKVhbHljVmHsStZIh1bIfzY5t3zJJ453IuMseOGbfH45bsMFvDt6qqnE0xSwkd9yZr5uyzCcmanG2ZfvP87GVxzJ9pzDLz1cDNha0Zn2OpKCNfJWdcfq7kuwleoYw9CtwBAQCCIAABAIIgAAEAgiAAAQCCIAABAIIojSy4sW7YlMcMO8ktc8ra8CuPnLMADeu31jGzNK9z/UxT3TNbFpgyjp85kJvxFtnmbcvgcswCzBzNbbDnnB2WzyZwPj7T19gBvrOhs8n+oAKdH3dAAIAgCEAAgCAIQACAIAhAAIAgCEAAgCBKIwvORT6zWxwzbZy6mfrqiOqSIeRQq871M2372rt8umVlZQYGTnnfyFLHzbVrq3EMx26etqzGWFm5Yat5DOfPNGx37XrrnKXoUu8wxHfWxiVrzjWTsEAz1fKJOyAAQBAEIABAEAQgAEAQBCAAQBDjLwnBlY9GU7ZGWy6JArahA5Rj8dKQzza2Y+kal/G9PVj3UOrFNXnEmJzgel154JyAYmHa30dyy7t/wEPJIR8JAeMwqcAVd0AAgCAIQACAIAhAAIAgCEAAgCAIQACAIEojC85HpppNiEyWUsueGeuGgZbPdM4O81B2xVxCR4qGBi37WzLBXLMdDeKVlcbtGUMpIms2nmUezvP2cX5CZKoV+HfTx3kYS9wBAQCCIAABAIIgAAEAgiAAAQCCIAABAIJwCkBf/epXFYvFRrwWLFiQfX9gYEAtLS2aNm2aJk2apFWrVqm7u9v7pHNEUe6rCMTKyoyv0e7rjem4FtqxdZijr2PoMkZ0Ysj4ss0lGh42vmJl5TkvxWLGl23szMCA8WU8fpaxbcfbNu+CuYZs68nnK4DoxAnjq2DOwwc43wGdddZZOnz4cPb1zDPPZN9bu3atHnvsMW3ZskXt7e06dOiQrrzySq8TBgCUBudf/8rKylRfX5+zPZVK6cEHH9TDDz+sSy+9VJK0adMmLVy4ULt27dJFF11kHC+dTiv9vn+H0Nvb6zolAEARcr4D2r9/v2bNmqUPf/jDuvrqq3Xw4EFJUkdHh4aGhtTc3Jzdd8GCBZozZ4527txpHa+trU3JZDL7amxsPI1lAACKjVMAWrp0qR566CE98cQT2rhxow4cOKCPf/zj6uvrU1dXlyoqKlRbWzviz9TV1amrq8s6Zmtrq1KpVPbV2dl5WgsBABQXp7+CW7FiRfa/Fy1apKVLl2ru3Ln6yU9+oqqqqtOaQCKRUCKROK0/CwAoXqNKo6qtrdUZZ5yhV199VZ/4xCc0ODionp6eEXdB3d3dxmdGcKt9lc8ul1au3SJDdKjMZ/dYh3GsHU4tn+lam2usO6K6rsdLRlWI6y2fY+MPGtW/Azp69Khee+01NTQ0aMmSJSovL9f27duz7+/bt08HDx5UU1PTqCcKACgtTndAf/d3f6fLL79cc+fO1aFDh3T77bdrwoQJ+uxnP6tkMqlrrrlG69at09SpU1VTU6M1a9aoqanJmgEHABi/nALQ73//e332s5/VW2+9pRkzZujiiy/Wrl27NGPGDEnS3XffrXg8rlWrVimdTmv58uW699578zJxAEBxi0VRYf2lZm9vr5LJpC7RSpXFzL1UMEZK7BmQlYexnZ+Z+JhLHr+6QfrK5PN6C3Etj2MnoiHt0KNKpVKqqamx7kctOABAEKXREbVIWX/LDJHxZhKiW2Q+f8MM8duu62c6zCVWXmHe1ZQxd7K5mHZ1vdNxPIama9/b3eJo9z3Z/uPlzsjlbxZGsXbugAAAQRCAAABBEIAAAEEQgAAAQZCE4JPjA0qnh66+Glz5SFHNJ9d1ujws9nUMDeM4J474KC00NOj2mS7ynIDidLwK6fq08ZEmX0jrHKPP5A4IABAEAQgAEAQBCAAQBAEIABAEAQgAEARZcD7lM3MoRIkaR15KC/la51gfr3wXuxzjYqRW+cxStO0botCpTRF8D4sJd0AAgCAIQACAIAhAAIAgCEAAgCAIQACAIMZfFlw+6y0VUi2nAHOxZruFqJPl4TOdsvp8HVcf9ery2abcV2v0sf5OjMMMsxEK9PxwBwQACIIABAAIggAEAAiCAAQACIIABAAIYvxlwQWo5eRcI81H1pjL/vnOmPMxjq9j6/KRto61Y31+ZF6ndY2+MtUcxCZMME8ln9c4Tl2BHlvugAAAQRCAAABBEIAAAEEQgAAAQRCAAABBFG4WXCyWmylToJkcWZYsI+eMrLHOYPPVtTMAH9luRbvOAury6ZTt5ukzMcbycC65AwIABEEAAgAEQQACAARBAAIABFG4SQhRJGkMH1S6PGAL8WA1xGeOk0Z6xfBA3FiKx1YqyCbEOovg2Dop1qSKfCYrjQJ3QACAIAhAAIAgCEAAgCAIQACAIAhAAIAgCjcLbqyFyAZxyUwptQw7G9fPLPTsRU+8lBzC6BVrZqiP7xVZcACAUkEAAgAEQQACAARBAAIABOEcgF5//XV97nOf07Rp01RVVaWPfexj2rNnT/b9KIp02223qaGhQVVVVWpubtb+/fu9ThoAUPycAtA777yjZcuWqby8XI8//rj27t2rb33rW5oyZUp2n7vuukv33HOP7rvvPu3evVvV1dVavny5BgYG3Gb2XkO69798MI3ra2xXUWR+meRz3i7zKDQucy/mdWL8CnHdjtHnOaVhf+Mb31BjY6M2bdqU3TZv3rzsf0dRpPXr1+srX/mKVq5cKUn64Q9/qLq6Oj3yyCP6zGc+42naAIBi53QH9NOf/lTnn3++Pv3pT2vmzJk699xz9cADD2TfP3DggLq6utTc3JzdlkwmtXTpUu3cudM4ZjqdVm9v74gXAKD0OQWg3/3ud9q4caPmz5+vbdu26frrr9cXv/hF/eAHP5AkdXV1SZLq6upG/Lm6urrsex/U1tamZDKZfTU2Np7OOgAARcYpAGUyGZ133nm64447dO655+raa6/VF77wBd13332nPYHW1lalUqnsq7Oz87THAgAUD6cA1NDQoDPPPHPEtoULF+rgwYOSpPr6eklSd3f3iH26u7uz731QIpFQTU3NiBcAoPQ5BaBly5Zp3759I7b99re/1dy5cyW9m5BQX1+v7du3Z9/v7e3V7t271dTU5DazfGVhFGsmVD7nbcuwy2fmXSFlI6Kwldq1UmrrGQWnLLi1a9fqj//4j3XHHXfoL//yL/Xss8/q/vvv1/333y9JisViuummm/S1r31N8+fP17x583Trrbdq1qxZuuKKK/IxfwBAkXIKQBdccIG2bt2q1tZW/eM//qPmzZun9evX6+qrr87u86UvfUn9/f269tpr1dPTo4svvlhPPPGEKisrvU8eAFC8YlFUWH//1Nvbq2QyqUu0UmWx8tDTGR9cb/99XDJF3BoBY6zUrpVSW4/BiWhIO/SoUqnUSZ/rUwsOABBE6TakGwe/ZXhTag32bErtmii19diwnpLFHRAAIAgCEAAgCAIQACAIAhAAIAgCEAAgiNLNgiPTJH9CZF/5+MxiviZM6y/m9QDiDggAEAgBCAAQBAEIABAEAQgAEETBJSG8Vxv1hIYknrEWqBAlYMZJ2RkrkhBQPE5oSNL//3luU3ABqK+vT5L0jP4j8ExgFeLn3nj/WTve14+i1NfXp2QyaX2/4NoxZDIZHTp0SJMnT1ZfX58aGxvV2dlZ0q26e3t7WWeJGA9rlFhnqfG9ziiK1NfXp1mzZiketz/pKbg7oHg8rtmzZ0t6t8OqJNXU1JT0yX8P6ywd42GNEussNT7XebI7n/eQhAAACIIABAAIoqADUCKR0O23365EIhF6KnnFOkvHeFijxDpLTah1FlwSAgBgfCjoOyAAQOkiAAEAgiAAAQCCIAABAIIgAAEAgijoALRhwwb90R/9kSorK7V06VI9++yzoac0Kk8//bQuv/xyzZo1S7FYTI888siI96Mo0m233aaGhgZVVVWpublZ+/fvDzPZ09TW1qYLLrhAkydP1syZM3XFFVdo3759I/YZGBhQS0uLpk2bpkmTJmnVqlXq7u4ONOPTs3HjRi1atCj7L8ebmpr0+OOPZ98vhTV+0J133qlYLKabbropu60U1vnVr35VsVhsxGvBggXZ90thje95/fXX9bnPfU7Tpk1TVVWVPvaxj2nPnj3Z98f6Z1DBBqAf//jHWrdunW6//XY9//zzWrx4sZYvX64jR46Entpp6+/v1+LFi7Vhwwbj+3fddZfuuece3Xfffdq9e7eqq6u1fPlyDQwMjPFMT197e7taWlq0a9cuPfnkkxoaGtInP/lJ9ff3Z/dZu3atHnvsMW3ZskXt7e06dOiQrrzyyoCzdjd79mzdeeed6ujo0J49e3TppZdq5cqVeuWVVySVxhrf77nnntP3v/99LVq0aMT2UlnnWWedpcOHD2dfzzzzTPa9UlnjO++8o2XLlqm8vFyPP/649u7dq29961uaMmVKdp8x/xkUFagLL7wwamlpyf7/8PBwNGvWrKitrS3grPyRFG3dujX7/5lMJqqvr4+++c1vZrf19PREiUQi+td//dcAM/TjyJEjkaSovb09iqJ311ReXh5t2bIlu8+vf/3rSFK0c+fOUNP0YsqUKdE//dM/ldwa+/r6ovnz50dPPvlk9Kd/+qfRjTfeGEVR6ZzL22+/PVq8eLHxvVJZYxRF0Ze//OXo4osvtr4f4mdQQd4BDQ4OqqOjQ83Nzdlt8Xhczc3N2rlzZ8CZ5c+BAwfU1dU1Ys3JZFJLly4t6jWnUilJ0tSpUyVJHR0dGhoaGrHOBQsWaM6cOUW7zuHhYW3evFn9/f1qamoquTW2tLToU5/61Ij1SKV1Lvfv369Zs2bpwx/+sK6++modPHhQUmmt8ac//anOP/98ffrTn9bMmTN17rnn6oEHHsi+H+JnUEEGoDfffFPDw8Oqq6sbsb2urk5dXV2BZpVf762rlNacyWR00003admyZTr77LMlvbvOiooK1dbWjti3GNf50ksvadKkSUokErruuuu0detWnXnmmSW1xs2bN+v5559XW1tbznulss6lS5fqoYce0hNPPKGNGzfqwIED+vjHP66+vr6SWaMk/e53v9PGjRs1f/58bdu2Tddff72++MUv6gc/+IGkMD+DCq4dA0pHS0uLXn755RF/n15KPvrRj+rFF19UKpXSv/3bv2n16tVqb28PPS1vOjs7deONN+rJJ59UZWVl6OnkzYoVK7L/vWjRIi1dulRz587VT37yE1VVVQWcmV+ZTEbnn3++7rjjDknSueeeq5dffln33XefVq9eHWROBXkHNH36dE2YMCEn06S7u1v19fWBZpVf762rVNZ8ww036Gc/+5l+/vOfZ/s7Se+uc3BwUD09PSP2L8Z1VlRU6CMf+YiWLFmitrY2LV68WN/5zndKZo0dHR06cuSIzjvvPJWVlamsrEzt7e265557VFZWprq6upJY5wfV1tbqjDPO0Kuvvloy51KSGhoadOaZZ47YtnDhwuxfN4b4GVSQAaiiokJLlizR9u3bs9symYy2b9+upqamgDPLn3nz5qm+vn7Emnt7e7V79+6iWnMURbrhhhu0detWPfXUU5o3b96I95csWaLy8vIR69y3b58OHjxYVOs0yWQySqfTJbPGyy67TC+99JJefPHF7Ov888/X1Vdfnf3vUljnBx09elSvvfaaGhoaSuZcStKyZcty/knEb3/7W82dO1dSoJ9BeUlt8GDz5s1RIpGIHnrooWjv3r3RtddeG9XW1kZdXV2hp3ba+vr6ohdeeCF64YUXIknRt7/97eiFF16I/vu//zuKoii68847o9ra2ujRRx+NfvWrX0UrV66M5s2bFx0/fjzwzE/d9ddfHyWTyWjHjh3R4cOHs69jx45l97nuuuuiOXPmRE899VS0Z8+eqKmpKWpqago4a3e33HJL1N7eHh04cCD61a9+Fd1yyy1RLBaL/vM//zOKotJYo8n7s+CiqDTWefPNN0c7duyIDhw4EP3iF7+Impubo+nTp0dHjhyJoqg01hhFUfTss89GZWVl0de//vVo//790Y9+9KNo4sSJ0b/8y79k9xnrn0EFG4CiKIq++93vRnPmzIkqKiqiCy+8MNq1a1foKY3Kz3/+80hSzmv16tVRFL2bBnnrrbdGdXV1USKRiC677LJo3759YSftyLQ+SdGmTZuy+xw/fjz627/922jKlCnRxIkToz//8z+PDh8+HG7Sp+Fv/uZvorlz50YVFRXRjBkzossuuywbfKKoNNZo8sEAVArrvOqqq6KGhoaooqIi+tCHPhRdddVV0auvvpp9vxTW+J7HHnssOvvss6NEIhEtWLAguv/++0e8P9Y/g+gHBAAIoiCfAQEASh8BCAAQBAEIABAEAQgAEAQBCAAQBAEIABAEAQgAEAQBCAAQBAEIABAEAQgAEAQBCAAQxP8DeHPi6CMJLN0AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(image_lens)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}