{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b7d08f1d",
   "metadata": {},
   "source": [
    "# Inverting the Lens Equation\n",
    "\n",
    "The lens equation $\\vec{\\beta} = \\vec{\\theta} - \\vec{\\alpha}(\\vec{\\theta})$ allows us to find a point in the source plane given a point in the image plane. However, sometimes we know a point in the source plane and would like to see where it ends up in the image plane. This is not easy to do since a point in the source plane may map to multiple locations in the image plane. There is no closed form function to invert the lens equation, in large part because the deflection angle $\\vec{\\alpha}$ depends on the position in the image plane $\\vec{\\theta}$. To invert the lens equation, we will need to rely on optimization and a little luck to find all the images for a given source plane point. Below we will demonstrate how this is done in caustic!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "4027aaf9",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "from functools import partial \n",
    "\n",
    "import torch\n",
    "from torch.nn.functional import avg_pool2d\n",
    "import matplotlib.pyplot as plt\n",
    "from ipywidgets import interact\n",
    "from astropy.io import fits\n",
    "import numpy as np\n",
    "from time import process_time as time\n",
    "\n",
    "import caustic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "2118e1c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# initialization stuff for an SIE lens\n",
    "\n",
    "cosmology = caustic.FlatLambdaCDM(name = \"cosmo\")\n",
    "cosmology.to(dtype=torch.float32)\n",
    "n_pix = 100\n",
    "res = 0.05\n",
    "upsample_factor = 2\n",
    "fov = res * n_pix\n",
    "thx, thy = caustic.get_meshgrid(res/upsample_factor, upsample_factor*n_pix, upsample_factor*n_pix, dtype=torch.float32)\n",
    "z_l = torch.tensor(0.5, dtype=torch.float32)\n",
    "z_s = torch.tensor(1.5, dtype=torch.float32)\n",
    "lens = caustic.SIE(\n",
    "    cosmology = cosmology, \n",
    "    name = \"sie\",\n",
    "    z_l = z_l,\n",
    "    x0 = torch.tensor(0.),\n",
    "    y0 = torch.tensor(0.),\n",
    "    q = torch.tensor(0.4),\n",
    "    phi = torch.tensor(np.pi/5),\n",
    "    b = torch.tensor(1.),\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "98e46aa1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Point in the source plane\n",
    "sp_x = torch.tensor(0.2)\n",
    "sp_y = torch.tensor(0.2)\n",
    "\n",
    "# Points in image plane\n",
    "x, y = lens.forward_raytrace(sp_x, sp_y, z_s)\n",
    "\n",
    "# Raytrace to check\n",
    "bx, by = lens.raytrace(x, y, z_s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "eb73147c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAGdCAYAAAAvwBgXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABRLElEQVR4nO3dd3iN9/8/8OfJRJYVIxJERe1Qo1QRpahSquhQQtWn39YKWqWqqqqhRs2qllpFFaFUUTP2CjFDJRKJ2CtLZJzcvz9evyROBYnknPfJOc/Hdd3XyblzJ+eVwf3Me+o0TdNAREREpICN6gKIiIjIejGIEBERkTIMIkRERKQMgwgREREpwyBCREREyjCIEBERkTIMIkRERKQMgwgREREpY6e6gCfJyMjAlStX4OLiAp1Op7ocIiIiygVN05CQkAAPDw/Y2Dy5zcOsg8iVK1fg5eWlugwiIiJ6BjExMfD09HziNWYdRFxcXADIF+Lq6qq4GiIiIsqN+Ph4eHl5Zd3Hn8Ssg0hmd4yrqyuDCBERUSGTm2EVHKxKREREyjCIEBERkTIMIkRERKSMWY8RISKibHq9HmlpaarLIAIA2Nvbw9bWNt+fh0GEiKgQSExMxOXLl6FpmupSiADIQFRPT084Ozvn6/MwiBARmTm9Xo/Lly+jWLFicHd35wKPpJymabh58yYuX74MHx+ffLWMMIgQEZm5tLQ0aJoGd3d3FC1aVHU5RAAAd3d3REVFIS0tLV9BhINViYgKCbaEkDkpqN9HBhEiIiJShkGEiIiMQtM0/O9//0PJkiWh0+kQGhqquqRnotPpsG7dOtVlWCwGESIiMorNmzdj0aJF+Ouvv3D16lXUrl1bdUmFSlRUVKEOcLnFwapERGQUERERKF++PF566aVn/hyapkGv18POzri3q9TUVDg4OJjd5zLH1ytobBEhIrISej2waxewYoU86vXGe60+ffpg0KBBiI6Ohk6nQ+XKlQEAKSkpGDx4MMqUKYMiRYrg5ZdfxpEjR7I+bteuXdDpdNi0aRMaNGgAR0dHbNy4Eba2tjh69CgAICMjAyVLlkSTJk2yPu63336Dl5dX1vPPP/8c1apVQ7FixVClShWMGTPGYDG4r7/+GvXq1cP8+fPh7e2NIkWKAAAuXLiAFi1aoEiRIqhZsya2bt361K/Vz88PAwcOREBAAEqXLo127doBAKZNm4Y6derAyckJXl5e+OSTT5CYmAgASEpKgqurK1avXm3wudatWwcnJyckJCTA29sbAFC/fn3odDr4+fllfW+7dOmCCRMmwMPDA88//zwAYOnSpWjYsCFcXFxQrlw5vPfee7hx44bB5z9z5gw6duwIV1dXuLi4oHnz5oiIiMh6//z581GjRg0UKVIE1atXx48//vjUrz+/2CJCRGQFgoKAIUOAy5ezz3l6AjNmAF27FvzrzZgxA8899xx+/vlnHDlyJGt654gRI7BmzRosXrwYlSpVwvfff4927dohPDwcJUuWzPr4kSNHYsqUKahSpQpKlCiBevXqYdeuXWjYsCFOnToFnU6H48ePIzExEc7OzggODkbLli2zPt7FxQWLFi2Ch4cHTp06hf79+8PFxQUjRozIuiY8PBxr1qxBUFAQbG1tkZGRga5du6Js2bI4dOgQ4uLiEBAQkKuvd/Hixfj444+xb9++rHM2NjaYOXMmvL29cfHiRXzyyScYMWIEfvzxRzg5OeGdd97BwoUL0a1bt6yPyXzu4uKCw4cPo3Hjxti2bRtq1apl0Oqxfft2uLq6GgSltLQ0jB8/Hs8//zxu3LiBYcOGoU+fPvj7778BALGxsWjRogX8/PywY8cOuLq6Yt++fUhPTwcALFu2DF999RVmz56N+vXr4/jx4+jfvz+cnJzg7++fq+/DM9HMWFxcnAZAi4uLU10KEZEyycnJ2tmzZ7Xk5ORn+vg1azRNp9M0wPDQ6eRYs6aAC/7/fvjhB61SpUpZzxMTEzV7e3tt2bJlWedSU1M1Dw8P7fvvv9c0TdN27typAdDWrVtn8LmGDRumvf7665qmadr06dO1t99+W/P19dU2bdqkaZqmVa1aVfv5558fW8vkyZO1Bg0aZD0fO3asZm9vr924cSPr3JYtWzQ7OzstNjY269ymTZs0ANratWsf+7lbtmyp1a9f/wnfCbFq1SqtVKlSWc8PHTqk2draaleuXNE0TdOuX7+u2dnZabt27dI0TdMiIyM1ANrx48cNPo+/v79WtmxZLSUl5Ymvd+TIEQ2AlpCQoGmapo0aNUrz9vbWUlNTc7z+ueee05YvX25wbvz48VrTpk1zvP5Jv5d5uX+za4aIyILp9dISktPK8JnnAgKM202TKSIiAmlpaWjWrFnWOXt7ezRu3BhhYWEG1zZs2NDgecuWLbF3717o9XoEBwfDz88Pfn5+2LVrF65cuYLw8PCsrgsAWLlyJZo1a4Zy5crB2dkZX375JaKjow0+Z6VKleDu7p71PCwsDF5eXvDw8Mg617Rp01x9bQ0aNHjk3LZt29C6dWtUqFABLi4u6NWrF27fvo379+8DABo3boxatWph8eLFAKR7qVKlSmjRosVTX69OnTqPjAsJCQlBp06dULFiRbi4uGS1EGV+3aGhoWjevDns7e0f+XxJSUmIiIhAv3794OzsnHV8++23Bl03xsAgQkRkwfbsMeyO+S9NA2Ji5Dpz4uTkZPC8RYsWSEhIwLFjx7B7926DIBIcHAwPDw/4+PgAAA4cOICePXuiQ4cO+Ouvv3D8+HGMHj0aqampT3yNgqw3KioKHTt2RN26dbFmzRqEhIRgzpw5AGBQx4cffohFixYBkG6Zvn375mqhsP++XlJSEtq1awdXV1csW7YMR44cwdq1aw1e70mr8maOXfnll18QGhqadZw+fRoHDx58aj35wSBCRGTBrl4t2Ovy47nnnoODg4PBOIq0tDQcOXIENWvWfOLHFi9eHHXr1sXs2bNhb2+P6tWro0WLFjh+/Dj++usvg/Eh+/fvR6VKlTB69Gg0bNgQPj4+uHTp0lPrq1GjBmJiYnD1oW/Gs96EQ0JCkJGRgalTp6JJkyaoVq0arly58sh177//Pi5duoSZM2fi7NmzBmMxMls89Llorjp37hxu376NiRMnonnz5qhevfojA1Xr1q2LPXv25LiDc9myZeHh4YGLFy+iatWqBkfmoFljYRAhIrJg5csX7HX54eTkhI8//hifffYZNm/ejLNnz6J///64f/8++vXr99SP9/Pzw7Jly7JCR8mSJVGjRg2sXLnSIIj4+PggOjoav//+OyIiIjBz5sys1oEnadOmDapVqwZ/f3+cOHECe/bswejRo5/pa61atSrS0tIwa9YsXLx4EUuXLsVPP/30yHUlSpRA165d8dlnn6Ft27bw9PTMel+ZMmVQtGhRbN68GdevX0dcXNxjX69ixYpwcHDIer3169dj/PjxBtcMHDgQ8fHxeOedd3D06FFcuHABS5cuxfnz5wEA48aNQ2BgIGbOnIl///0Xp06dwsKFCzFt2rRn+h7kFoMIEZEFa95cZsc8rrVfpwO8vOQ6U5g4cSLeeust9OrVCy+88ALCw8OxZcsWlChR4qkf27JlS+j1eoOxIH5+fo+ce+ONNzB06FAMHDgQ9erVw/79+zFmzJinfn4bGxusXbsWycnJaNy4MT788ENMmDDhWb5M+Pr6Ytq0aZg0aRJq166NZcuWITAwMMdr+/Xrh9TUVHzwwQcG5+3s7DBz5kzMmzcPHh4e6Ny582Nfz93dHYsWLcKqVatQs2ZNTJw4EVOmTDG4plSpUtixYwcSExPRsmVLNGjQAL/88kvWmJEPP/wQ8+fPx8KFC1GnTh20bNkSixYtMnqLiE7TchrCZB7i4+Ph5uaGuLg4uLq6qi6HiEiJBw8eIDIy0mC9i7wICgIyZ4g+/D9+ZjhZvdo4U3gpd5YuXYqhQ4fiypUrhWphsif9Xubl/s0WESIiC9e1q4SNChUMz3t6MoSodP/+fURERGDixIn46KOPClUIKUgMIkREVqBrVyAqCti5E1i+XB4jIxlCVPr+++9RvXp1lCtXDqNGjVJdjjLsmiEiMnP57ZohMgZ2zRAREVGhxyBCREREyjCIEBERkTIMIkRERKQMgwgREREpY9QgEhgYiEaNGsHFxQVlypRBly5dspaSJSIiIjJqEAkODsaAAQNw8OBBbN26FWlpaWjbti2SkpKM+bJEREQm9fXXX6NevXqqyyiU7Iz5yTdv3mzwfNGiRShTpgxCQkLQokULY740ERERFQJGDSL/lblzYMmSJXN8f0pKClJSUrKex8fHm6QuIiKLFhcHJCTImu7/dfky4OICuLmZvq4CptfrodPpYGPD4Y+Ficl+WhkZGQgICECzZs1Qu3btHK8JDAyEm5tb1uHl5WWq8oiILFNcHNC+PdCyJRATY/i+mBg53769XFfAVq9ejTp16qBo0aIoVaoU2rRpk9U1n5GRgW+++Qaenp5wdHREvXr1DFrRd+3aBZ1Oh3v37mWdCw0NhU6nQ1RUFABpZS9evDjWr1+PmjVrwtHREdHR0UhJScHnn38OLy8vODo6omrVqliwYEHW5zl9+jRee+01ODs7o2zZsujVqxdu3br12K8j83XWrVsHHx8fFClSBO3atUPMf7+fDzly5AheffVVlC5dGm5ubmjZsiWOHTtmcI1Op8P8+fPx5ptvolixYvDx8cH69esNrslrrYWRyYLIgAEDcPr0afz++++PvWbUqFGIi4vLOp70QyYiolxISABu3AAuXgT8/LLDSEyMPL94Ud6fkFCgL3v16lW8++67+OCDDxAWFoZdu3aha9euyNxVZMaMGZg6dSqmTJmCkydPol27dnjjjTdw4cKFPL3O/fv3MWnSJMyfPx9nzpxBmTJl0Lt3b6xYsQIzZ85EWFgY5s2bB2dnZwDAvXv38Morr6B+/fo4evQoNm/ejOvXr6NHjx5PfZ0JEyZgyZIl2LdvH+7du4d33nnnsdcnJCTA398fe/fuxcGDB+Hj44MOHTog4T/f53HjxqFHjx44efIkOnTogJ49e+LOnTv5qrXQ0UxgwIABmqenp3bx4sU8fVxcXJwGQIuLizNSZURE5i85OVk7e/aslpyc/GyfIDpa06pU0TRAHvftM3weHV2wBWuaFhISogHQoqKicny/h4eHNmHCBINzjRo10j755BNN0zRt586dGgDt7t27We8/fvy4BkCLjIzUNE3TFi5cqAHQQkNDs645f/68BkDbunVrjq87fvx4rW3btgbnYmJiNADa+fPnc/yYzNc5ePBg1rmwsDANgHbo0CFN0zRt7Nixmq+vb44fr2maptfrNRcXF23Dhg1Z5wBoX375ZdbzxMREDYC2adOmZ67VlJ70e5mX+7dRW0Q0TcPAgQOxdu1a7NixA97e3sZ8OSIiyomXF7BrF1ClirSANGsmj1WqyHkjdIP7+vqidevWqFOnDrp3745ffvkFd+/eBSDj/65cuYJmzZoZfEyzZs0QFhaWp9dxcHBA3bp1s56HhobC1tYWLVu2zPH6EydOYOfOnXB2ds46qlevDgCIiIh47OvY2dmhUaNGWc+rV6+O4sWLP7be69evo3///vDx8YGbmxtcXV2RmJiI6Ohog+sert3JyQmurq64ceNGvmotbIw6WHXAgAFYvnw5/vzzT7i4uODatWsAADc3NxQtWtSYL01ERA/z8gKWLpUQkmnpUqOEEACwtbXF1q1bsX//fvzzzz+YNWsWRo8ejUOHDqFUqVJP/fjMAafaQxvEp6WlPXJd0aJFodPpDJ4/SWJiIjp16oRJkyY98r7y5cs/ta7c8vf3x+3btzFjxgxUqlQJjo6OaNq0KVJTUw2us7e3N3iu0+mQkZFh0lpVM2qLyNy5cxEXFwc/Pz+UL18+61i5cqUxX5aIiP4rJgbo1cvwXK9ejw5gLUA6nQ7NmjXDuHHjcPz4cTg4OGDt2rVwdXWFh4cH9u3bZ3D9vn37ULNmTQCAu7s7ABlrkik0NPSpr1mnTh1kZGQgODg4x/e/8MILOHPmDCpXroyqVasaHE5OTo/9vOnp6Th69GjW8/Pnz+PevXuoUaNGjtfv27cPgwcPRocOHVCrVi04OjrmeZDps9Za2Bi9ayano0+fPsZ8WSIietjDA1OrVAH27cvupnl4AGsBOnToEL777jscPXoU0dHRCAoKws2bN7Nu3J999hkmTZqElStX4vz58xg5ciRCQ0MxZMgQAEDVqlXh5eWFr7/+GhcuXMDGjRsxderUp75u5cqV4e/vjw8++ADr1q1DZGQkdu3ahT/++AOAtNTfuXMH7777Lo4cOYKIiAhs2bIFffv2hV6vf+zntbe3x6BBg3Do0CGEhISgT58+aNKkCRo3bpzj9T4+Pli6dCnCwsJw6NAh9OzZM889Ac9aa2HDydZERJbs8mXDELJrF/DSS4ZjRvz85LoC5Orqit27d6NDhw6oVq0avvzyS0ydOhWvvfYaAGDw4MEYNmwYhg8fjjp16mDz5s1Yv349fHx8AMiNf8WKFTh37hzq1q2LSZMm4dtvv83Va8+dOxfdunXDJ598gurVq6N///5Z04YzW2L0ej3atm2LOnXqICAgAMWLF3/i+iPFihXD559/jvfeew/NmjWDs7PzE1v3FyxYgLt37+KFF15Ar169MHjwYJQpUya337581VrY6LSHO+DMTHx8PNzc3BAXFwdXV1fV5RARKfHgwQNERkbC29sbRYoUydsHZ64jcuPGowNTM1tKypQBNm+2iEXNjGHRokUICAgwWNOEnvx7mZf7t0lXViUiIhNzc5OQkdPKql5eQHCwxaysSoUTgwgRkaVzc3t80Mhp2XciE7KcTiYiIiIj6NOnD7tljIhBhIiIiJRhECEiIiJlGESIiIhIGQYRIiIiUoZBhIiIiJRhECEiIiJlGESIiIhIGQYRIiIyCj8/PwQEBKgug8wcV1YlIrIS+gw99kTvwdWEqyjvUh7NKzaHrY2t0po0TYNer4edHW9H1ootIkREViAoLAiVZ1RGq8Wt8F7Qe2i1uBUqz6iMoLAgo7xenz59EBwcjBkzZkCn00Gn0yEqKgq7du2CTqfDpk2b0KBBAzg6OmLv3r3o06cPunTpYvA5AgIC4Ofnl/U8IyMDgYGB8Pb2RtGiReHr64vVq1c/sY7KlStj/PjxePfdd+Hk5IQKFSpgzpw5BtdER0ejc+fOcHZ2hqurK3r06IHr169nvf/EiRNo1aoVXFxc4OrqigYNGuDo0aP5/h6RYBAhIrJwQWFB6PZHN1yOv2xwPjY+Ft3+6GaUMDJjxgw0bdoU/fv3x9WrV3H16lV4PbTz78iRIzFx4kSEhYWhbt26ufqcgYGBWLJkCX766SecOXMGQ4cOxfvvv4/g4OAnftzkyZPh6+uL48ePY+TIkRgyZAi2bt0KQMJN586dcefOHQQHB2Pr1q24ePEi3n777ayP79mzJzw9PXHkyBGEhIRg5MiRsLe3f4bvCuWEbWFERBZMn6HHkM1DoEF75H0aNOigQ8DmAHR+vnOBdtO4ubnBwcEBxYoVQ7ly5R55/zfffINXX301158vJSUF3333HbZt24amTZsCAKpUqYK9e/di3rx5aNmy5WM/tlmzZhg5ciQAoFq1ati3bx9++OEHvPrqq9i+fTtOnTqFyMjIrKC0ZMkS1KpVC0eOHEGjRo0QHR2Nzz77DNWrVwcA+Pj45Lpuejq2iBARWbA90XseaQl5mAYNMfEx2BO9x4RVAQ0bNszT9eHh4bh//z5effVVODs7Zx1LlixBRETEEz82M7g8/DwsLAwAEBYWBi8vL4PWmpo1a6J48eJZ1wwbNgwffvgh2rRpg4kTJz719Shv2CJCRGTBriZcLdDrCoqTk5PBcxsbG2iaYatNWlpa1tuJiYkAgI0bN6JChQoG1zk6OhqpSvH111/jvffew8aNG7Fp0yaMHTsWv//+O958802jvq61YBAhIrJg5V3KF+h1eeHg4AC9Xp+ra93d3XH69GmDc6GhoVljMWrWrAlHR0dER0c/sRsmJwcPHnzkeY0aNQAANWrUQExMDGJiYrJaRc6ePYt79+6hZs2aWR9TrVo1VKtWDUOHDsW7776LhQsXMogUEHbNEBFZsOYVm8PT1RM66HJ8vw46eLl6oXnF5gX+2pUrV8ahQ4cQFRWFW7duISMj47HXvvLKKzh69CiWLFmCCxcuYOzYsQbBxMXFBZ9++imGDh2KxYsXIyIiAseOHcOsWbOwePHiJ9axb98+fP/99/j3338xZ84crFq1CkOGDAEAtGnTBnXq1EHPnj1x7NgxHD58GL1790bLli3RsGFDJCcnY+DAgdi1axcuXbqEffv24ciRI1lBhvKPQYSIyILZ2thiRvsZAPBIGMl8Pr39dKOsJ/Lpp5/C1tYWNWvWhLu7O6Kjox97bbt27TBmzBiMGDECjRo1QkJCAnr37m1wzfjx4zFmzBgEBgaiRo0aaN++PTZu3Ahvb+8n1jF8+HAcPXoU9evXx7fffotp06ahXbt2AACdToc///wTJUqUQIsWLdCmTRtUqVIFK1euBADY2tri9u3b6N27N6pVq4YePXrgtddew7hx4/L53aFMOu2/nXJmJD4+Hm5uboiLi4Orq6vqcoiIlHjw4AEiIyPh7e2NIkWKPNPnCAoLwpDNQwwGrnq5emF6++noWqNrQZVqdipXroyAgACu8GoET/q9zMv9m2NEiIisQNcaXdH5+c5mt7IqEYMIEZGVsLWxhV9lP9VlEBlgECEiIosVFRWlugR6Cg5WJSIiImUYRIiIiEgZBhEiokLCjCc5khUqqN9HBhEiIjNnayszW1JTUxVXQpQt8/cx8/fzWXGwKhGRmbOzs0OxYsVw8+ZN2Nvbw8aGf0OSWhkZGbh58yaKFSsGO7v8RQkGESIiM6fT6VC+fHlERkbi0qVLqsshAiAbFVasWBE6Xc7bB+QWgwgRUSHg4OAAHx8fds+Q2XBwcCiQ1jkGESKiQsLGxuaZl3gnMlfsaCQiIiJlGESIiIhIGQYRIiIiUoZBhIiIiJRhECEiIiJlGESIiIhIGQYRIiIiUoZBhIiIiJRhECEiIiJlGESIiIhIGQYRIiIiUoZBhIiIiJRhECEiIiJlGESIiIhIGTvVBRCR5dDrgT17gKtXgfLlgebNAVtb1VURkTljECGiAhEUBAwZAly+nH3O0xOYMQPo2lVdXURk3tg1Q0T5FhQEdOtmGEIAIDZWzgcFqamLiMwfgwgR5YteLy0hmvbo+zLPBQTIdURE/8UgQkT5smfPoy0hD9M0ICZGriMi+i8GESLKl6tXC/Y6IrIuDCJElC/lyxfsdURkXRhEiChfmjeX2TE6Xc7v1+kALy+5jojovxhEiChfbG1lii7waBjJfD59OtcTIaKcMYgQUb517QqsXg1UqGB43tNTznMdESJ6HC5oRkQFomtXoHNnrqxKRHnDIEJEBcbWFvDzU10FERUm7JohIiIiZdgiQkRkIdLSgLt3gXv3gPh4ICEBuH8fSE4GUlJkddv0dBlEbGsL2NkBRYrI4eICuLoCJUoAZcrIOSJTYBAhIjJzmgZcvw5ERgLR0XJcvgxcuSLjca5fB27ckPBRUFxdZbBxxYqAtzdQvbocdesC5coV3OsQGTWI7N69G5MnT0ZISAiuXr2KtWvXokuXLsZ8SSKiQisuDjh3To7z54F//5UjIkJaNnLL1VUOFxfAyUlaNxwdpQUkc/BwRoa0oKSkyOdOSJAgc+eOnI+PB86eleO/PDyABg1kMLKfH/DCCxyUTM/OqEEkKSkJvr6++OCDD9CV8/eIiABIF8m//wKhoXKcPAmcOSN78jyOjY20UFSqJK0Unp4yXbp8eaBsWelOKV0aKF48f6FA0yQQXbsm9URHA+HhEozCwqTuK1fk2LBBPqZ0aaBDB5k11aEDu3Uob3SaltOemUZ4IZ0uzy0i8fHxcHNzQ1xcHFxdXY1XHBGRkWRkABcuAEeOAEePAiEhwPHjQFJSztd7eGR3g/j4ANWqyWOlSoCDg2lrz0lSkoSngweB4GBg924JLpnc3IAePYABAwBfX2VlkmJ5uX+bVRBJSUlBSkpK1vP4+Hh4eXkxiBBRoZGUBBw6BOzdCxw4IG/fvfvodcWKyY26Xj0Zd1GnDlCrlrRoFCZpacD+/cD69cCqVYatOq++CowYAbRu/fgtAMgy5SWImNVg1cDAQIwbN051GUREuXbnjizitnu3hI9jx2RmysOKFAHq1wcaNQIaNpQxFdWrW8a4Cnt7oGVLOSZPllaSefMklGzdKkerVrINQJ06qqslc8QWESKiPEhKktCxbRuwfbuM7/jv/6KenjKQs1kzoEkTafGwt1dTrypRUcC0acDPP8uAWFtbYNgw4JtvOIbEGhTaFhFHR0c4OjqqLoOIKIteL2M7tmyR8HHwoHRHPKxGDWkRePllCR+VKysp1axUrgzMnAkMHy7HmjXSYrJpExAUJONeiAAzCyJERObg1i1g82Zg40bgn3+k++VhlSrJ+IfWraXboWxZNXUWBpUqycaH69cD//sfcPq0dFEFBQGvvKK6OjIHRg0iiYmJCA8Pz3oeGRmJ0NBQlCxZEhUrVjTmSxMR5ZqmASdOAH/9JeHj0CHD7hY3N6BNGwkfbdoAVapw8GVevfGGBJDu3YF9+4D27SWctG+vujJSzahjRHbt2oVWrVo9ct7f3x+LFi166sdz+i4RGUt6ugwyXbtWboiXLhm+v25d4PXX5XjxRVkMjPLvwQPg/felq6ZIEenyatFCdVVU0Mxy+u6zYBAhooKUnCyzONaskcW4Hp5WW7SotHh07CiLclWooK5OS5eaCrz1lrRAlSolrVH8fluWQjtYlYiooCUny1/dq1ZJy0diYvb7SpeW4PHmm9LlUqyYujqtiYMD8McfMrD3+HGgd28JiDbcD94qMYgQkcVJSZHBpitXSsvHw+HDywvo2lWOZs0sYy2PwqhoUWDFCllTZccOYNEi4IMPVFdFKrBrhogsgl4vi2n9/rvM0ni426ViRekK6NFDxntwoKn5mDIF+OwzWdr+wgW2SlkKds0QkdU4eRJYskT+ur5yJfu8hwfw9ttyNG7M8GGuBg0C5syRBdCWLgU++kh1RWRqDCJEVOhcuwYsXy4B5MSJ7PMlSkjLx7vvygJj7HYxf46OEkaGDwd++olBxBqxa4aICoXUVBls+uuvMvg0I0POOzgAnToBvXrJmhRcnLnwuXMHKFdOVqz991+uumoJ2DVDRBbj1CkJH7/9JiueZmraVMLH228DJUuqq4/yr2RJ2Ztnxw5ZAp5BxLowiBCR2UlMlEGnP/8MHDmSfb58eaBPH6BvX96sLE3r1hJEDh4EBg9WXQ2ZEoMIEZmNY8ckfCxfDiQkyDl7e1kevG9foF07rnBqqerVk8eTJ5WWQQrwnzQRKZWcLOt9/PijYetH1aqySZq/P1CmjLr6yDSqVpXH6Gi1dZDpMYgQkRLh4cDcucDChdlrftjby6yX/v1lV1tOubUemWEzIUH2oylSRG09ZDoMIkRkMhkZspT3jBkyKDFTpUrA//0f0K8f4O6urj5Sp2jR7LdTUxlErAmDCBEZXWKitHzMmiWrZwLS2tG+PTBggDxyzQ/rZr4LSZCxMYgQkdFcuiThY/58IC5Ozrm6ysDTgQOzxwUQZf5+6HSAk5PaWsi0GESIqMAdOABMmwYEBWUvPFatmkzL9PcHnJ3V1kfm58YNeSxRgq1j1oZBhIgKREaG7HQ7eTKwb1/2+TZtgIAA4LXXuM07Pd7Fi/Lo7a22DjI9BhEiypcHD4DFi4GpU7PHfzg4AO+/DwwdCtSurbY+KhzOnpXHatXU1kGmxyBCRM/k3j1g3jxg+nTZhA4AiheX2S+DB8sqqES5dfSoPDZooLYOMj0GESLKkxs3gB9+kK3bM1c/9fIChg2T6bcuLmrro8JH04C9e+Xtxo3V1kKmxyBCRLkSEyPdLz//LKuhAkCtWsCIEcC778piZETP4tQpCbjFigEvvqi6GjI1BhEieqKLF4HAQBkHkpYm5xo1AkaPBjp14gBUyr/16+WxVSsZX0TWhUGEiHIUHg5MmAAsXQro9XLOzw/44guZCcPl1wsPvR7Yswe4elXG7jRvbl5TZFetkseuXdXWQWowiBCRgYgIYPx44LffsgNI27bAV18BzZqprY3yLigIGDIEuHw5+5ynpyyzbw43/pMn5bC3Bzp3Vl0NqcBGVSICIF0w/foBzz8v3TB6PdChA3DwILBlC0NIYRQUBHTrZhhCACA2Vs4HBamp62GLFsljx45AqVJKSyFFGESIrNzly8BHH0kA+fVXCSDt20sA2biRgwcLK71eWkJy2sMl81xAQHarlwr378seRADwwQfq6iC1GESIrNSNG3Kjeu45mQmTni5dMPv3y864DCCF2549j7aEPEzTZCbUnj2mq+m/fvtN1qOpXFlW3iXrxDEiRFYmPl6m4U6dCiQlybmWLYFvvwVeflltbVRwrl4t2OsKWno68P338vaQIeY1eJZMi0GEyEo8eCCLkAUGArdvy7mGDeV569acBWNpcruyraoVcH//XQZGlyoF9O+vpgYyD+yaIbJwer0MPq1WDfj0Uwkh1asDq1cDhw9zKq6lat5cZsc87mer08mKuM2bm7YuQFpDvvlG3h4+HHByMn0NZD4YRIgslKbJWI969YA+fWQ8gKcnsGCBrGT51lsMIJbM1lam6AKP/pwzn0+frqZL5NdfZYPEUqWAQYNM//pkXhhEiCxQaKgMPO3QATh9WjajmzQJ+PdfmZ1gx05Zq9C1q7R8VahgeN7TU86rWEckIUHWpAHk0dnZ9DWQeeF/R0QW5OpV4MsvZUqkpsly2YMHy2qoJUqoro5U6NpVFgozl5VVv/kGuH4dqFpVdmomYhAhsgDJyTILZuLE7Jkw77wDfPcd4O2ttjb6j4QEoHt3GSn81Vcm2VzF1laW51ctLEy6gwB55L4yBLBrhqhQ0zRgxQpZjGzMGAkhTZoABw7IeYYQM3L9OtClC+DqKkvVTpgAREWprspkMjKkBSQ9XTZLfP111RWRuWAQISqkQkKkif2992QgasWKEj7275cwQmZA04AjR4CePYFy5YA//8x+X+fOMpXJSixYAOzeDRQrBsyapboaMifsmiEqZG7cAEaPlv/YNU3+Yx81SqZBFi2qujoCANy5AyxbBsyfLzu6/VfmFCYrERUFDBsmb3/7LVCpktJyyMwwiBAVEunpwNy50gUTFyfnevaU2TD/nRVBCqSkAP/8I81SQUHy/L9Wr5Z501YkI0NmaiUmysq9gwerrojMDYMIUSGwfz/wySfAiRPyvF49YPZs7oirXHo6sGOHLBO6dq1snJKTDh1km1l3d1NWZxZ++AHYuVNa7n79lUu506MYRIjM2K1bwIgR2TuUlighYxz/9z/+h65MUhKwdSvw118y5uPWrez3lS8vI4d37ZLn9vbSZDVkCGBjfUPyQkJk6jgggcTHR209ZJ4YRIjMUEaG/PX4+ecy3ACQ5u2JE03wR/UffwCOjjKYkkR0NLBxI7Bhg7SAPNztUrq0TMd9+23gyhWgXz85X6UKsHKlTNO1QvHx8i1JTZXJQtxPhh6HQYTIzJw8KdMcDxyQ576+MjakaVMTvHhystw9APnzdedO6xyAcueOfO3bt0vwOH/e8P3e3jIHtVMnWaBDp5PBO4GB8v62bWWsSMmSJi/dHGia5LGICJnN9euv3E6AHo9BhMhM3L8vq05OmSIb1Tk7y/NBg0y4JPvD024uXJAw8sknclSpYqIiFIiNleR34IAEkNBQuZtmsrGRJNixo4SPmjWz76xJScC770prCSA7CwYGWvU6+j/8IONy7e2lUYir+tKT6DTt4X9t5iU+Ph5ubm6Ii4uDq6ur6nKIjOaff6QVJDJSnr/5JjBzpqIZnrGxj76wTge0by/LtXbqVLjvLImJsuvfoUPZ4SMm5tHratYEWrcGXnkFaNky56/5xg0JJ0eOAEWKyHTdnj2N/zWYse3bgXbtJEzPng0MGKC6IlIhL/dvBhEihW7flvUVliyR556ewJw5wBtvqK0L3boBa9ZIv1D58sDmzdnvs7OT7ohWreQG3aiRea7VnZ4OXLoku/6dOJF9REQ8eq2NDVC3rqwE9/LLEj7Kl3/y5w8Pl3AWESFdMBs2AC+9ZJyvpZCIjJRfh9u3gd69ZaIQu2SsU17u39bbdkikkKYBq1YBAwcCN2/Kf9aDBsliTy4uqquDtHqsWSN/5W/aJN00y5bJudOngW3b5ACkO8fXF6hTR27mdesCzz0nN3JjzhTRNLnjxcbKINHISKkz84iMBNLScv5YDw/ghReku6VpU7l75mUb2KNHZUruzZsyXmTTJpktY8Xi4qRx6PZtGZ87bx5DCOUOW0SITOzqVRlysW6dPK9VS1r0zWpZ9vBwGR/i4CCDVx6eK3z+vOyVEhwsa3Y/PH31Yfb20sRTqRJQpoy0GmQerq7yuTMPe3uZKpSWln2kpsr4i7t3ZX2Oe/fk7Tt3JHhcuZLzomEPK1JEAoKvr+FRuvSzf2/27wdee02mhbzwgsymKVfu2T+fBUhLk+y6ZYtkvMOHrXOMM2Vj1wyRGdI0aVQYPFjup3Z2slT7F1+YYc9GerqsQJWWJt0bFSvmfJ2mSTA5eTL7OH1aprvq9aap1d1d7npeXhKeHj48PQu2VSY4WHZrS0qSbqkNG8ykCUsdTZPxTT//LL8ye/ZIPiPrxq4ZIjNz7Rrw0UfA+vXy/IUXZJGyunXV1vVYdnbyV35MjBT/uCCi0wHVq8vRo0f2+fR0abG4dElCye3b0pKReSQkZLd6ZB42NtIy8vDh7AwULy4DRYsXz37bw0PCR7lysuaJKWzbJoN3kpOBNm1kMbNixUzz2mZs4kQJITqdzFhmCKG8YhAhMrKVK6Ur5s4dubeOHSurpdrbq67sKTJvssnJef9YOzsJL48LMAXlwAHg/Hnoe/fCnug9uJpwFeVdyqN5xeawXbJUumUKYgGWXbuk7+HBAxkbkjl+xsotXGi4cqryQdZUKDGIEBnJnTsSQFaulOf16wOLF8uYzkIhcx2M9HS1dTzOgQPASy9BA/DZ3wH4oVZc1ruGnnHD1FVx0AEypiM/YeTIkewQ0rGjLJBhqlYYM7Z+ffZqqZ9/LqvYEz0L69v8gMgEtmyRwLFypYzzHDMGOHiwEIUQQAZjAoCbm9o6Huf8eWgAdACmroqDf4ic9g9BVgjR/v91z+z0aZmim5goU3pXrWIIgSw226OHDAPy989eUJboWTCIEBWg5GSZhtu+vQyReP55+cP9m2/McEDqk2ha9iY3ZhpE9L17YXh3t6wwsnADsGS1PGaGkE+7F4e+d69ne4HISFmq/c4d4MUXZUwIu2Nw4IB0waSkyMJ78+dzmi7lD4MIUQE5cULWT5g9W54PGAAcOyZLVBQ6t27JzBBA0fKuT7cneg9+qBWHvp2QFUZ6nc4OIX07AdNq3cOe6D15/+R378pYkKtXpRnr77/zts6IhTp+XGYuJyXJeN0VK6x6JXsqIAwiRPmUkSED9Ro3Bs6eBcqWlfvW7NmFeFLFv//KY8WKhvvPmJGrCVcBAIsbAL/VNnzfb7Xl/MPX5Vpqqqwse+6chLDNm61287qHhYZK+IiLk8Vn161jLxUVDAYRonzI3Gpk2DC5f3XqJNuYvPaa6sry6fhxeaxRQ20dT1DeRZZg9w8B3j9t+L73TyNrzEjmdbmSuSjGjh3SArJxo0wVtnInT0oIyeyl2rgRcHJSXRVZCgYRome0fbss0rlpkwwdmDNHhhG4u6uurADs2yePzZqpreMJmldsjqFn3AzGhCytDYMxI8POFEfzis1z/0lnzpQ5qTY2MtLYbBd6MZ1jx2Rbodu3pZtxyxZZGJeooDCIEOVRejrw5ZfAq6/KWl+1askMz08+sZBBe5oG7N0rb5txELFdstRgdkzfTkDvbjAYMzJl1T1ZTyQ39u0DPv1U3p46VcaIWLkjR2QD4jt3pOvxn3/MduwyFWIMIkR5cOWK/Mc8YYLcr//3P9lXo3btp39soXHiBHD5sowNMasNcP7j+eezQsjw7m5ZY0IWN5DZMplhJFeb0V2/LvNR09OBt9/mohiQbYRat5Ytfl56Cdi6VRa2JSpoHO9MlEvbtgHvvScbrjo7A7/8ArzzjuqqjODPP+WxbVvzHm3btCmwfz90589jcu9eeOM/K6vqOuRyZVW9Hnj3XUmZNWpwPipksPVbb8kabq1aya+ElW+pQ0bEIEL0FBkZ0gIydqy0gvj6yrpWPj6qKzMCTZMvDgA6d1ZbS240bQo0bQpbAH6V/Qzf16dP7j7H1KnAzp2SLoOCrH6a7ooVQO/e0jjUsaP8OnD5FDImds0QPcGdO/Kf8VdfyT36ww9lQSeLDCGADAo4c0buPG++qboa4ztxQgb8AMCMGbJ5nxX78UegZ08JIe+9J7mMIYSMjUGE6DFCQ2WBssxZMQsXSneMmS6rUTDmz5fHbt0sf0DAgwfA++/LLsCdOwN9+6quSBlNk7A9YIC8/cknwNKlhWBjRrII7JohysHy5dL6kZwMeHvLX4b16qmuysju3ZN2eQDo109pKSYxfrzsJVOmTPY+9lYoLQ34+GNgwQJ5PnasHFb67SAFTNIiMmfOHFSuXBlFihTBiy++iMOHD5viZYnyLD1dFifr2VNCyGuvAUePWkEIAYC5c2Vzt9q1gZYtVVdjXGfPApMny9tz50oYsUKJidIYtGCBLJ0ydy7w9dcMIWRaRg8iK1euxLBhwzB27FgcO3YMvr6+aNeuHW7cuGHslybKkzt3JHj88IM8Hz0a2LDBSlb3fvBAxkgAwIgRln0n0jRpAkhLk6VwrWEsTA5iY4EWLaTrsWhRWbL9//5PdVVkjYweRKZNm4b+/fujb9++qFmzJn766ScUK1YMv/76q7FfmijXzp6VBZu2bZMZq6tXA99+C9jaqq7MRObPl7U0vLwsdE7yQ5Ytk0UyihUDZs2y7ND1GCdOyFLtx4/LSsA7d0omI1LBqEEkNTUVISEhaNOmTfYL2tigTZs2OHDgwCPXp6SkID4+3uAgMrYtW2QWaEQEULmyzIp56y3VVZlQfDzwzTfy9qhRlj1C8cEDaeoCZLZMpUpq61Fg/XpZMDc2VpZNOXRIQgmRKkYNIrdu3YJer0fZsmUNzpctWxbXrl175PrAwEC4ubllHV5eXsYsjwizZ8tK3vHxsqPo4cNWuL3I1KmySlu1ajJC15LNmQNERwMVKgABAaqrMSlNk2ExXboASUmyauq+fTIYm0gls5q+O2rUKMTFxWUdMTExqksiC6XXyyregwbJgmV9+ki3jEVsWJcXsbESRADgu+8suzXk3j1ZmQ6QFiCLnodtKDlZFikbMSJ7iMymTUCJEqorIzLy9N3SpUvD1tYW169fNzh//fp1lCtX7pHrHR0d4ejoaMySiJCUJCt6b9ggzwMDgc8/t8qhAsDQofINadoU6NpVdTXGNWcOcPcuULMm4O+vuhqTiY2VH+3hwzLmafp0YOBA1VURZTNqi4iDgwMaNGiA7du3Z53LyMjA9u3b0fRp+z8QGcGNG8Arr0gIcXQE/vgDGDnSSkPIpk2yfretrczbtORvQnJy9qygL76wmlHIe/cCDRpICClZUnbPZQghc2P0Bc2GDRsGf39/NGzYEI0bN8b06dORlJSEvla8iiGpER4OtGsHXLwIlColg/Zeekl1VYrcvy/LaALSR+Xrq7YeY1u0SMbBVKoku+taOE2TbBkQILOU69SR6blVqqiujOhRRg8ib7/9Nm7evImvvvoK165dQ7169bB58+ZHBrASGdOxY7JGyI0bMjhv06bc7Q5vsT77DIiMBDw9gXHjVFdjXJom3TKArFZnZ9kLSt+/L+uBLF0qz3v0AH79FXByUlsX0ePoNE3TVBfxOPHx8XBzc0NcXBxcXV1Vl0OF1M6dwBtvyCqS9epJCMlhiJL12LRJpgoB0lb/6qtq6zG2Q4eAJk1kcOqVKxa9h86FC0D37rJOiK0tMGmSZC9L7nUj85SX+7dl/2lAVu/PP6UlPiUF8POT5mk3N9VVKXT7NvDBB/L24MGWH0IAq9nILyhI9u2Lj5fZXytXAq1aqa6K6OnMavouUUFatkwWJktJkbUTNm2y8hCSkSF3qmvXZLv7iRNVV2R8aWmyTC6QHcAsTGqqTH56663s9XCOH2cIocKDQYQs0rx5QK9esl5I794yOaRIEdVVKTZxYvZ0oeXLrWMdjd27Zf0Qd3egeXPV1RS4qCj5sqZPl+fDhwM7dsh6bUSFBYMIWZyZM2WwnqbJxJCFCy1+fOLTbdsGjBkjb8+ZA9Svr7YeU/nzT3ns1MnipuyuWSNjng4floXJ1q8Hpkyx7DXpyDIxiJBFmTJFZqMCsorkrFmyvblVi46WFdwyMoB+/eSwFplrGL3+uto6ClBysqyM2q0bEBcn43CPH+emdVR4Wft/0WRBpkyRWamA7Gc2cSJnCyA+Xm7Ct25JK8isWaorMp27d2VbZUAGTliAU6eAhg2Bn36S5yNHSu+TFe7dRxbE2husyUL88EN2CPn6a2DsWKXlmIf0dJkydPq0zFdet846xoVkOnRIHn18gDJl1NaST5omGzR+9pkMvi5bFliyBGjbVnVlRPnHIEKF3ty5slYCIMMgGEIgd67Bg4HNmyV8bNgAVKyouirTymwNqVdPaRn5de2aTHbavFmev/66LFBWyLMVURZ2zVChtmQJ8Mkn8vbIkZa/SGiufftt9v4xy5dLe761+fdfeaxWTW0d+bBunSzPvnmzzPqaOVMyJUMIWRK2iFChtXat/KUIAIMGyS72Vj8mBJDN3b76St7+4QdZRMUaXb4sj5UrKy3jWcTHyz4xCxfKc19fWRenVi2lZREZBVtEqFDasQN45x2ZCNKnj6yjwBACuXMFBMjbX3+dPYXIGt27J48lSigtI6927gTq1pUfpU4HfP65DHdhCCFLxRYRKnRCQ4E335QVJd96C/jlF07RBQD88Qfw4Yfy9tCh2a0i1io+Xh4LyT5VSUnAqFHZE5u8vYHFiy1yHTYiAwwiVKhERQHt28s9pmVLaa62+sXKAOC33wB//6y1QvTfT8WeYB2uXgXKl5ebmYWt5/V0mSt7paerrSMX9uyRbsaICHn+0UfA5MmAi4vaujLp9VKjVf8+kdHwv3AqNO7dk01jr1+Xpus//5TVyq3er79KS4imAf36Iaj9PAzx1mUNkQAAT08ZOtK1awG+bmqq7GYbGyvHrVvyQ7p7Vx7v3ZMVt1JT5UhLk8f0dLmLOToCDg5yFCkiXSglS2Yf5cvLAhmVKsl81bw2e2VOVb5/vwC/6IKVmAiMHi2tIJomP6f584F27VRXli0oSHr4jP77RFaLQYQKhbQ0WUkyLAzw8AA2brTyDewyzZ2bPW3o448R9MpsdOthA00zvCw2Vr5/q1fn4eahafIn8IUL2Ud4OHDxonzCmzcL9Et5IkdHWQ+kTh1JoXXrAo0ayR4yj1OunDw+fAc1I9u3S36MipLnH34oi/KZ0+91UJD83hTI7xPRYzCIkNnTNGDgQPmP28lJQoinp+qqFNM0GYz6zTfyPCAA+snTMMRb98hNI/NynU7GsXbunEOzekICcPIkcOJE9nH6tAxceBIHB9lhzcND5pSWKCFH8eJyuLlJiLC3l2vt7aUvLSMju6UkJUXWLb93D7hzR45bt+Rud+mSPKakSD2nTwMrVmS/fs2a0kf3yivSjPBwX0bVqvJ4/nxevrNGd/cu8Omn0pAFyPIuv/xifouT6fXSEvJMv09EecAgQmbvxx+Bn3+W//hWrCj061PlX1qa7OqXeScbMwYYNw57gnVP/ONf04CYGGDPznT4FQ8FDhyQ4/Dh7MEJ/2VjI9NffXyyj+eeA7y8JHyUKmX86UppadKqcfashKVTp2TEcliYnDt7VlqGHB0ljHTvLn+uZ66dsnevcevLJU2TFoRBg6R7EZDGrIkTzWcsyMP27HlyY1LW79MewM/PZGWRBWIQIbMWHJw9G/X777mxFxITgR49gE2bJCT8+KOMbIT0ouTG1Q79gLQlj76jQgVZsKJu3ezHqlWlJUMle3uZQuLtbbh53a1bchcMDgb+/lu6jtavlyPzT3VAgktEhAQoRaKjJXRs3CjPq1eXsSDNmikr6aly/fuUy+uIHodBhMxWbKz8cZueLpvHDh+uuiLFLl2SxclCQ2Ug5sqVBsmsfPncfZryaZek+6RJE6BpU3msXx8oXdooZRtN6dIyj/vNN2XhttOnpclh4UL5Uz2zxQgApk0D5swxeYlpabIa6tix0stlby9TdL/4wvwHWuf69ymX1xE9jk7TcuoBNA/x8fFwc3NDXFwcXAvJWgBUMNLSpLl3/375w/zAAaBYMdVVKbRrl6SyW7dkgOaGDcCLL2a///p16Df8jcqDOiH2QUloOaxVqEMGPEveR+SeWNjWqGa5K8ClpwN//QVMmgQcPJh9/swZGVNiIgcOSA/ayZPy/OWXgXnzTFpCvuj10isXG5vzOBGdTsZqRUZyjAg9Ki/3by4DRWZpzBgJIW5uwJo1VhxCMrddbdNGQsgLLwAhIRJCzp4FJkyQFo3y5WHb/wPMePA/ABI6HqbTAdDZYPovzrCt+bzlhhBABsN26SK/QKtXZ5/fsMEkL3/7NtC/P/DSSxJCSpWSxpng4MITQgAJFzNmyNv//XXJfD59OkMIFQDNjMXFxWkAtLi4ONWlkAlt2qRpcgfWtNWrVVejUEKCpvXqlf3N6NlT006c0LTx4zWtdu3s85lHgwaa9vXX2prJEZqnZ4bBu7y8NG3NGtVfkCJXr2razz9r2u3bRn0ZvV7TfvlF00qVyv6+9+2raTdvGvVljW7NGk3z9NT4+0R5kpf7N7tmyKzcuCFLRdy4AQwYII0BVun4ceDtt2UAJiALelWoABw7ln2Nvb3M+XzjDaBjR5nF8v9xJUzTCgmRwaiHD8vz2rVlIs/LL6utq6Dw94nyKi/3bwYRMhuaJosjrVsnG3wdPSoLblqVjAwZFHPmTM7vt7WVbpq335buh0K2oZuluXVLBp7Ony+/vy4usrzLoEHZK8wTWaO83L85a4bMxvLlEkLs7WUPGasLIRcvPn6KadOmQO/essvfk1YTJZNIS5MWj6+/lgXKAKBnT5li/lDDFBHlAoMImYWrV2X1VEA2jfX1VVuPSWVkyPreCxcani9ZEujTB+jXr3CNcrRwW7bI5sZhYfLc11e6EC2lG4bI1BhEyCwMGSIrfDdsCIwcqboaEzpwQKZXPMzVVdb87tJF/WJilCUsTNay2bRJnpcuLZOW+vXjeAmi/GAQIeW2bAFWrZL/zH/5RWZfWo3x4w2fHzxouD4IKXfrlnTB/PSTDNq0t5fWu6++ku10iCh/rOm/fDJDKSkysA8ABg+2wn1kRo6UnfwCAsx7vW8r9OABMGuWtHrExcm5zp2ByZNlyx0iKhgMIqTUjBkyQ7VcOVkG2+q0aCEHmY2MDFk9f9QoWVUfkIA8bRrQqpXS0ogsEldWJWVu3ZK/NgHZgdTNTW09RNu2AY0aAe+9JyGkQgUZQ3z0KEMIkbGwRYSUmTgRiI+XvzZ79VJdDVmzkBBZD+Sff+S5iwswYgQwbJgVby9AZCIMIqTEtWuygz0AfPed7GhPZGrnz8u+RqtWyXN7e+Djj4Evv+RyLUSmwiBCSkyZAiQny35t7durroasTUwM8M030u2i18smbu+/D4wbB3h7q66OyLowiJDJ3b0r26ED8teoJW8ES+bl2jVpgZs3D0hNlXOdOslYpTp11NZGZK0YRMjkFiwAEhNlS5XXXlNdDVmDGzdk+fUff5SWOADw85MA8t/15IjItBhEyKQyMrLHhgwezNYQMq6bN2XdjzlzgPv35dyLL0oAad1abW1EJBhEyKSCg4HISJmq++67qqshS3Xjhqz7MXs2kJQk5xo1kjEg7dszABOZEwYRMqklS+Tx7bc5LZIK3rVrMhB67tzsFpAGDSSAdOjAAEJkjhhEyGRSU4G1a+Xtnj3V1kKWJSpKumAWLJBtAwDZQPGrr4COHRlAiMwZgwiZzM6dsmdH2bLcVoUKxvnzQGAgsGwZkJ4u55o0kQDCLhiiwoFBhEwmc9XKjh25bTrlz5EjwKRJQFAQoGlyrk0bYPRooGVLBhCiwoRBhExm5055bNNGbR1UOGkasGWLdMHs2JF9/o03JIA0bqyuNiJ6dgwiZBL37wMnT8rbL7+sthYqXFJTgRUrgKlTgVOn5JydnWxM99lnQO3aausjovxhECGTOHVKltIuW1Z2NCV6mrt3gV9+AWbMAK5ckXPOzkD//sCQIUClSmrrI6KCwSBCJvHvv/JYowb77+nJzp+X8LF4cfYU3PLlZQG8jz4CSpRQWx8RFSwGETKJqCh5rFJFaRlkpjQN2LYNmD4d+Pvv7PN16wIBAdIN4+ioqjoiMiYGETKJ27flsUwZtXWQeUlIkEXuZs2SlhBAWsw6dgSGDpX9YNiCRmTZGETIJBIS5NHVVW0dZB7CwmTPocWLs383nJ2BPn2kC8bHR2l5RGRCDCJkUvzr1nqlpQHr1skGdMHB2eerVQMGDgT8/RlUiawRgwiZRGb/fuYGZGQ9Ll4E5s8Hfv0VuH5dztnYAJ06SQBp3ZoBlciaMYiQSWSODblxQ20dZBqpqcCGDcBPP8kg1ExlywIffiizX7y81NVHROaDQYRMInPNh/BwtXWQcZ05Iy0fS5cCN2/KOZ0OePVV4H//k1VQ7e3V1khE5oVBhEyiTh15DA0FMjKkaZ4sw717wMqVEkAOH84+X64c8MEHsgBZ5cqqqiMic8cgQiZRty7g5ATcuQOcOAHUr6+6IsqPtDTZxHDJEuDPP4GUFDlvZydTb/v1k91v7fg/DBE9Bf+bIJNwcJBBievXy8wJBpHCR9OAY8eA334Dli83HO9TqxbQty/QqxfXiiGivGEQIZN55x0JIgsWAGPG8K/lwuLffyV4rFiRvVQ/ALi7y4qn/v5AvXqc+UJEz4a3AjKZrl2B0qWB2Fi5qfXqpboiepxLl4BVq4DffwdCQrLPFyki02579wbatePAUyLKPwYRMhlHR2D4cGDUKGD0aOCtt4BixVRXRZkyw8fq1cChQ9nnbW1l1st77wFdugAuLspKJCILpNM0TVNdxOPEx8fDzc0NcXFxcOWSixYhORl4/nkgJgYYMACYPVt1Rdbt/Hlg7VogKAg4ciT7vE4HtGwJdO8OdOvGcR9ElDd5uX+zRYRMqmhRWWWzXTtZ6rt1a+DNN1VXZT0yMmSK7Z9/SgDJ3GgOkPDRogXQo4d0o5Urp65OIrIeDCJkcm3bAsOGAdOmAT17Aps3yw2QjCMxUVY33bgR+Osv4Nq17PfZ2wOtWkk3WefOsvIpEZEpMYiQEhMnAufOAX//Dbz+uiwH7uenuirLoGkyu+XvvyV87N4t635kcnEBXntNWqJeew1wc1NXKxERgwgpYW8vgyJffx3YuVNaSX78UfYhoby7dUu+j9u2AVu2yMDThz33nHyvX39dAp+Dg5IyiYgeYbQgMmHCBGzcuBGhoaFwcHDAvXv3jPVSVEgVLSp/sffuLaGkf39gxw4ZwFqypOrqzFtiIrB/P7B9u4SP48elJSSTg4N0d3XoIEe1alzng4jMk9GCSGpqKrp3746mTZtiwYIFxnoZKuSKFpV9SiZOlEXOVqwAtm4FAgNlpU5bW9UVmoe7d4F9+4A9e6Sr5ehRID3d8JratWXwb9u2MuPFyUlNrUREeWH06buLFi1CQEDAM7WIcPqudTl8WMLH2bPyvHZtYORImcVhTQtnaRpw4QJw8KC0euzbB5w+/eh1FSsCr7wCtGkjj+XLm75WIqKcFNrpuykpKUjJ3D0L8oWQ9WjcWHbnnTMHGDdObr7vvy9hZMAAWYm1QgXVVRYsTZPxHCEhchw9Ksfdu49eW60a0Lw58PLLMs6DO9oSkSUwqyASGBiIcePGqS6DFLK3BwICZP+SuXOBmTOBy5dlNdYvvpAbcJcussy4t7fiYvPo/n0gLEwC1smTErpOnABu33702iJFgAYNgCZNJHi89BIXFSMiy5SnrpmRI0di0qRJT7wmLCwM1atXz3qel66ZnFpEvLy82DVjxVJSZNzIggXA3r2G7/PxkWDy8sty037+efUb6aWlAdHRQESEdK9cuCCLhp07Jy0fOf1rs7cH6tQBXngBaNhQjjp1OLOFiAqvvHTN5CmI3Lx5E7dz+vPtIVWqVIHDQ/+DcowIFZSoKGDNGllzZO9eQK83fH/RokDNmkD16tKNUbmyjKPw8JDWBDe3Z585kpoq3SW3bwM3bwLXr8vCYLGxcsTESACJiXm0roeVLi0ho3Zt2bHW11fednR8trqIiMyR0caIuLu7w93dPV/FET2rypVl07zhw4F79ySMBAcDBw5IN0dSUvZYi5zY2QGurhJInJyk+8PBQWbm2NjI8ufp6dKqkZwsR2IikJAgb+eWoyNQpYqEIR8fOWrUkIDEfz5ERIaM1pAdHR2NO3fuIDo6Gnq9HqGhoQCAqlWrwtnZ2VgvS1aieHGgY0c5AGmFCA+XMRjnzkmXSHS0HNeuAfHxEjLu3JHjWeh08rplyshRtqwMnvXwkJaXihWBSpVk9oqNTUF9pUREls1o03f79OmDxYsXP3J+586d8MvlWt7smqGC8uCBrD4aFydHcrKcS0mRcRt6vbSM2NnJUbSotJg4O0srSmZLCtc1ISJ6OqONETE1BhEiIqLCJy/3bzYgExERkTIMIkRERKQMgwgREREpwyBCREREyjCIEBERkTIMIkRERKQMgwgREREpwyBCREREyjCIEBERkTIMIkRERKQMgwgREREpwyBCREREyjCIEBERkTIMIkRERKQMgwgREREpwyBCREREyjCIEBERkTIMIkRERKQMgwgREREpwyBCREREyjCIEBERkTIMIkRERKQMgwgREREpwyBCREREyjCIEBERkTIMIkRERKQMgwgREREpwyBCREREyjCIEBERkTIMIkRERKQMgwgREREpwyBCREREyjCIEBERkTIMIkRERKQMgwgREREpwyBCREREyjCIEBERkTIMIkRERKQMgwgREREpwyBCREREyjCIEBERkTIMIkRERKQMgwgREREpwyBCREREyjCIEBERkTIMIkRERKQMgwgREREpwyBCREREyjCIEBERkTIMIkRERKQMgwgREREpwyBCREREyjCIEBERkTIMIkRERKQMgwgREREpwyBCREREyjCIEBERkTIMIkRERKQMgwgREREpwyBCREREyjCIEBERkTIMIkRERKQMgwgREREpY7QgEhUVhX79+sHb2xtFixbFc889h7FjxyI1NdVYL0lERESFjJ2xPvG5c+eQkZGBefPmoWrVqjh9+jT69++PpKQkTJkyxVgvS0RERIWITtM0zVQvNnnyZMydOxcXL17M1fXx8fFwc3NDXFwcXF1djVwdERERFYS83L+N1iKSk7i4OJQsWfKx709JSUFKSkrW8/j4eFOURURERIqYbLBqeHg4Zs2ahY8++uix1wQGBsLNzS3r8PLyMlV5REREpECeg8jIkSOh0+meeJw7d87gY2JjY9G+fXt0794d/fv3f+znHjVqFOLi4rKOmJiYvH9FREREVGjkeYzIzZs3cfv27SdeU6VKFTg4OAAArly5Aj8/PzRp0gSLFi2CjU3usw/HiBARERU+Rh0j4u7uDnd391xdGxsbi1atWqFBgwZYuHBhnkIIERERWT6jDVaNjY2Fn58fKlWqhClTpuDmzZtZ7ytXrpyxXpaIiIgKEaMFka1btyI8PBzh4eHw9PQ0eJ8JZwwTERGRGTNaX0mfPn2gaVqOBxERERHAvWaIiIhIIQYRIiIiUoZBhIiIiJRhECEiIiJlGESIiIhIGQYRIiIiUoZBhIiIiJRhECEiIiJlGESIiIhIGQYRIiIiUoZBhIiIiJRhECEiIiJlGESIiIhIGQYRIiIiUoZBhIiIiJRhECEiIiJlGESIiIhIGQYRIiIiUoZBhIiIiJRhECEiIiJlGESIiIhIGQYRIiIiUoZBhIiIiJRhECEiIiJlGESIiIhIGQYRIiIiUoZBhIiIiJRhECEiIiJl7FQX8CSapgEA4uPjFVdCREREuZV53868jz+JWQeRhIQEAICXl5fiSoiIiCivEhIS4Obm9sRrdFpu4ooiGRkZuHLlClxcXKDT6VSX80zi4+Ph5eWFmJgYuLq6qi7H6vHnYT74szAv/HmYl8L+89A0DQkJCfDw8ICNzZNHgZh1i4iNjQ08PT1Vl1EgXF1dC+Uvk6Xiz8N88GdhXvjzMC+F+efxtJaQTBysSkRERMowiBAREZEyDCJG5ujoiLFjx8LR0VF1KQT+PMwJfxbmhT8P82JNPw+zHqxKRERElo0tIkRERKQMgwgREREpwyBCREREyjCIEBERkTIMIiYSFRWFfv36wdvbG0WLFsVzzz2HsWPHIjU1VXVpVmvChAl46aWXUKxYMRQvXlx1OVZnzpw5qFy5MooUKYIXX3wRhw8fVl2SVdq9ezc6deoEDw8P6HQ6rFu3TnVJViswMBCNGjWCi4sLypQpgy5duuD8+fOqyzI6BhETOXfuHDIyMjBv3jycOXMGP/zwA3766Sd88cUXqkuzWqmpqejevTs+/vhj1aVYnZUrV2LYsGEYO3Ysjh07Bl9fX7Rr1w43btxQXZrVSUpKgq+vL+bMmaO6FKsXHByMAQMG4ODBg9i6dSvS0tLQtm1bJCUlqS7NqDh9V6HJkydj7ty5uHjxoupSrNqiRYsQEBCAe/fuqS7Farz44oto1KgRZs+eDUD2lfLy8sKgQYMwcuRIxdVZL51Oh7Vr16JLly6qSyEAN2/eRJkyZRAcHIwWLVqoLsdo2CKiUFxcHEqWLKm6DCKTSk1NRUhICNq0aZN1zsbGBm3atMGBAwcUVkZkXuLi4gDA4u8TDCKKhIeHY9asWfjoo49Ul0JkUrdu3YJer0fZsmUNzpctWxbXrl1TVBWRecnIyEBAQACaNWuG2rVrqy7HqBhE8mnkyJHQ6XRPPM6dO2fwMbGxsWjfvj26d++O/v37K6rcMj3Lz4OIyNwMGDAAp0+fxu+//666FKOzU11AYTd8+HD06dPniddUqVIl6+0rV66gVatWeOmll/Dzzz8buTrrk9efB5le6dKlYWtri+vXrxucv379OsqVK6eoKiLzMXDgQPz111/YvXs3PD09VZdjdAwi+eTu7g53d/dcXRsbG4tWrVqhQYMGWLhwIWxs2CBV0PLy8yA1HBwc0KBBA2zfvj1rUGRGRga2b9+OgQMHqi2OSCFN0zBo0CCsXbsWu3btgre3t+qSTIJBxERiY2Ph5+eHSpUqYcqUKbh582bW+/hXoBrR0dG4c+cOoqOjodfrERoaCgCoWrUqnJ2d1RZn4YYNGwZ/f380bNgQjRs3xvTp05GUlIS+ffuqLs3qJCYmIjw8POt5ZGQkQkNDUbJkSVSsWFFhZdZnwIABWL58Of7880+4uLhkjZlyc3ND0aJFFVdnRBqZxMKFCzUAOR6khr+/f44/j507d6ouzSrMmjVLq1ixoubg4KA1btxYO3jwoOqSrNLOnTtz/Hfg7++vujSr87h7xMKFC1WXZlRcR4SIiIiU4SAFIiIiUoZBhIiIiJRhECEiIiJlGESIiIhIGQYRIiIiUoZBhIiIiJRhECEiIiJlGESIiIhIGQYRIiIiUoZBhIiIiJRhECEiIiJlGESIiIhImf8HAXqVbqhkXmAAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "A = lens.jacobian_lens_equation(thx, thy, z_s)\n",
    "detA = torch.linalg.det(A)\n",
    "CS = ax.contour(thx, thy, detA, levels = [0.], colors = \"b\", zorder = 1)\n",
    "# Get the path from the matplotlib contour plot of the critical line\n",
    "paths = CS.collections[0].get_paths()\n",
    "caustic_paths = []\n",
    "for path in paths:\n",
    "    # Collect the path into a descrete set of points\n",
    "    vertices = path.interpolated(5).vertices\n",
    "    x1 = torch.tensor(list(float(vs[0]) for vs in vertices))\n",
    "    x2 = torch.tensor(list(float(vs[1]) for vs in vertices))\n",
    "    # raytrace the points to the source plane\n",
    "    y1,y2 = lens.raytrace(x1, x2, z_s)\n",
    "\n",
    "    # Plot the caustic\n",
    "    ax.plot(y1,y2, color = \"r\", zorder = 1)\n",
    "ax.scatter(x, y, color = \"b\", label = \"forward raytrace\", zorder = 10)\n",
    "ax.scatter(bx, by, color = \"r\", marker = \"x\", label = \"source plane\", zorder = 9)\n",
    "ax.scatter([sp_x.item()], [sp_y.item()], color = \"g\", label = \"true pos\", zorder = 8)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0aa515f",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "PY39",
   "language": "python",
   "name": "py39"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}